- Author:
- Ting Yu <ting.yu@auckland.ac.nz>
- Date:
- 2014-05-05 11:11:02+12:00
- Desc:
- Change the unit from dimensionless to TODO.
- Permanent Source URI:
- http://models.cellml.org/workspace/194/rawfile/c560868b736d95fd544954a4a3d5875175c194dc/ikeda_julie_16Feb07.txt
METHOD RK4
STARTTIME = 0
STOPTIME=60*3 {minutes}
DT = 0.01
;*****************************************************************
; 14 May 2010 (note by S. Randall Thomas):
; This is a Berkeley Madonna implementation of the Ikeda model of acid-base regulation
; Reference: Ikeda, N., F. Marumo, M. Shirataka and T. Sato (1979).
; "A model of overall regulation of body fluids."
; Annals of Biomedical Engineering 7: 135-166.
;
; This implementation was done by Julie Fontecave-Jallon at Univ. Joseph Fourier (Grenoble)
; in the SAPHIR project, which was also seedEP1 in the VPH Network of Excellence (FP7).
;
; Please treat this code as a privileged communication, since it is not yet published material.
;
;*****************************************************************
FCOI=0
YGLI=0
QVIN=0
QIN=0.001
;QIN=IF time>5 AND time<=10 THEN 0.2 ELSE 0 ;Fig 10
;YGLI=IF time>5 AND time<=55 THEN 1000 ELSE 0 ;Fig 12
;FCOI=IF time>5 AND time<=35 THEN 0.05 ELSE 0 ;Fig 11
;*****************************************************************
; BLOCK 1 - CARDIOVASCULAR SYSTEM
; Input = VB
; Output = QCO, PAS, PVS, PVP, PAP
;*****************************************************************
TRSP=2
; constant parameters
RTOT=20
KR=0.3
RTOP=3
KL=0.2
DEN=1
QCO0=VB*DEN
QCO=QCO0+1
PAS=20+RTOT*QCO0
PVS=MAX(0,-10.33+QCO0/KR)
PAP=8+RTOP*QCO0
PVP=MAX(-16+QCO0/KL,0)
;*****************************************************************
; BLOCK 2 - RESPIRATION
; Input = STBC, QCO, VTW
; Output = DCLA, PCOA, XCO3, PHA
;*****************************************************************
; constant parameters
XHB=15
FO2I=0.21
PBA=760
PBL=PBA-47
VAL=3
MRCO=0.2318
MRO2=0.2591
init (FCOA)=0.0561
init(FO2A)=0.1473
init(UCOV)=0.6075
init(UO2V)=0.1515
; equations (3) to (8)
d/dt(FCOA)=(VI*(FCOI-FCOA)+863/(PBA-47)*QCO*(UCOV-UCOA))/VAL
d/dt(FO2A)=(VI*(FO2I-FO2A)+863/(PBA-47)*QCO*(UO2V-UO2A))/VAL
d/dt(UCOV)=(MRCO+QCO*(UCOA-UCOV))/VTW
d/dt(UO2V)=(-MRO2+QCO*(UO2A-UO2V))/VTW
UCOA=6.732*10^-4*PCOA+0.02226*XCO3
UO2A=3.168*10^-5*PO2A+UHBO
; equation (1)
;VI=VR*VI0
d/dt(VI)=(VR*VI0-VI)/TRSP
init(VI)=5
; equation (2) et F21
k1=IF PHA<=7.4 THEN 0.22 ELSE 0.0258
k6=IF PHA<=7.4 THEN -12.734 ELSE -5.003
k3=0.58
k4=3.496
k2=IF PCOA>40 THEN 1 ELSE 0.0396
k5=IF PCOA>40 THEN -32.08 ELSE 160.11
VR=k1*H+k2*(k3+k4/(PO2A-32))*(PCOA+k5)+k6
LIMIT VR >= 0
VI0=5
H=10^(9-PHA)
; F23
f=(1-exp(-PO2A*g))^2
g=0.0066815*PHA^3-0.10098*PHA^2+0.44921*PHA-0.454
UHBO=UHB*f
; and also...
DCLA=XCO3-STBC
UHB=XHB/75
PCOA=FCOA*(PBA-47)
PO2A=FO2A*(PBA-47)
; F24
;XCO3=STBC-(0.527*XHB+3.7)*(PHA-7.4)+0.375*(UHB-UHBO)/0.02226
; equation (9) and F22
PHA=6.1+log10(XCO3/(0.03*PCOA))
XCO0=STBC-(0.527*XHB+3.7)*(PHA-7.4)+0.375*(UHB-UHBO)/0.02226
d/dt(XCO3)= (XCO0-XCO3)/TRSP
init(XCO3)=24
;*****************************************************************
; BLOCK 3 - EXTRACELLULAR SPACE
; Input = PAS, PVS, QIC, QWU
; Output = PPCO, VB; VEC, XPP
;*****************************************************************
; constant parameters
;QIN=0.001
;QVIN=0
QMWP=0.0005
QLF0=0.002
VIF0=8.8
QIWL=0.0005
CFC=0.007
VRBC=1.8
; en plus
VB=VRBC+VP
HT=VRBC/VB
VEC=VP+VIF
d/dt(VIN)=QIN-(VIN/10)
init (VIN)=0.01
d/dt(VIF)=QCFR-QLF-QIC
init (VIF)=VIF0
PC=(CRAV*PVS+PAS)/(1+CRAV)
CRAV=5.93
d/dt(ZPP)=YPLF-YPLG-YPLV-YPLC
init(ZPP)=154
XPP=ZPP/VP
YPLC=QPLC*(XPP-XPIF)
XPIF=ZPIF/VIF
d/dt(ZPIF)=YPLC-YPG-YPLF
init(ZPIF)=176
YPLF=XPIF*QLF
YPLV=XPP*0.00047-0.0329
YPLG=0.00023*(XPP-ZPLG)
d/dt(ZPLG)=(XPP-ZPLG)/24
init(ZPLG)=70
YPG=0.0057*(XPIF-ZPG)
d/dt(ZPG)=(XPIF-ZPG)/150
init (ZPG)=20
; equation (10)
VP_plus=VIN/10+QVIN+QMWP+QLF
; equation (11)
VP_moins=QIWL+QWU+QCFR
; equation (12)
d/dt(VP)=VP_plus-VP_moins
init(VP)=2.2
; equation (13)
QCFR=CFC*Pf
; equation (14)
Pf=PC-PPCO-PIF+PICO
; F31
x=VIF/VIF0
PIF=IF (x<=0.9) THEN (-15) ELSE IF (x>0.9 AND x<=1) THEN (87*x-93.3) ELSE IF (x>1 AND x<=2) THEN (-6.3*(2-x)^10) ELSE (x-2)
; F32
QLF=QLF0*(24/(1+exp(-0.4977*PIF)))
; F33
QPLC=2.768*10^-6*PC^2
; F34
PPCO=0.4*XPP
; F35
PICO=0.25*XPIF
;*****************************************************************
; BLOCK 4 - INTRACELLULAR SPACE and ELECTROLYTES
; Input = GFR, PHA, VEC, YKU, YNU
; Output = QSMP, QIC, VTW, XKE, XNE, YGLU, YHI, YMNU, YNU
;*****************************************************************
; constant parameters
CSM=0.0003
YNIN=0.12
CKEI=0.001
YKIN=0.047
XGLO=108
YINS=0
;YGLI=0
;init YGLI=0
;flow YGLI= (pulse(1, 5,STOPTIME) + pulse(-1, 55,STOPTIME))*1
YMNI=0
YURI=0.15
CGL1=1
CGL2=1
CGL3=0.03
CHEI=5
CBFI=10^-9
; equation (15)
d/dt(ZNE)=YNIN-YNU+YHI
init(ZNE)=1540
; equation (16)
d/dt(ZKE)=YKIN-YKU-y4
init(ZKE)=49.5
; en plus
d/dt(ZKI)=y4
init(ZKI)=2800
y4=z4+CKEI*(2800*F41-ZKI)
; F41
F41=1+0.5*log10(XKE/(56.744-7.06*PHA))
; F42
YGLU=IF XGLE*GFR<0.65 THEN 0 ELSE XGLE*GFR-0.65
; en plus
XNE=ZNE/VEC
XKI=ZKI/VIC
XKE=ZKE/VEC
d/dt(ZHI)=YHI
YHI=CHEI*(0.4-PHA+PHI)
init(ZHI)=100
PHI=-log10(CBFI*ZHI)
d/dt(YINT)=1/1.50*(XGLE-XGLO/18-YINT)
init(YINT)=0
YGLS=CGL1*YINT+CGL2*YINS
z4=CGL3*YGLS
d/dt(ZGLE)=YGLI/180-YGLS-YGLU
init(ZGLE)=66
XGLE=ZGLE/VEC
OSMP=(XNE+XKE)*1.86+XGLE+XURE+XMNE+9.73
d/dt(ZMNE)=YMNI-YMNU
init(ZMNE)=0
XMNE=ZMNE/VEC
d/dt(ZURE)=YURI-YURU
init(ZURE)=77.5
XURE=ZURE/VTW
VTW=VEC+VIC
YMNU=1*GFR*XMNE
YURU=XURE*GFR*0.6
d/dt(VIC)=QIC
init(VIC)=20
QIC=CSM*((-XNE-XKE)-XGLE+(10.5+XKI))
;*****************************************************************
; BLOCK 5 - KIDNEY 2
; Input = DCLA,GFR,PCOA,VEC,XCO3,XKE,XNE,XPP,YHI,YKU,YNH4,YNU, STPG
; Output = STBC, YCO3, YORG, YPO4
;*****************************************************************
; constant parameters
YCAI=0.007
YCLI=0.1328
YMGI=0.008
YOGI=0.01
YPOI=0.025 ;mM/min
;YPOI=0.045; mEq/min
YSOI=0.02
; F50
F50=-PCOA/120+4/3
; F51
YCO3=IF XCO3*GFR*F50<=2 THEN 0 ELSE IF (XCO3*GFR*F50>2 AND XCO3*GFR*F50<=4) THEN 0.1638*( XCO3*GFR*F50-2)^2.61 ELSE XCO3*GFR*F50-3
; F52
YCA=IF XCAE*GFR<0.493 THEN 0 ELSE XCAE*GFR-0.493
; F53
YMG=IF XMGE*GFR<0.292 THEN 0 ELSE XMGE*GFR-0.292
; F54
YSO4=IF XSO4*GFR<0.08 THEN 0 ELSE XSO4*GFR-0.08
; F55
YPO4=IF XPO4*GFR<=0.11 THEN 5/22*XPO4*GFR ELSE XPO4*GFR-0.085
;YPO4=IF XPO4*GFR<=0.198 THEN 1/8*XPO4*GFR ELSE XPO4*GFR-0.175
; F56
YORG=IF XOGE*GFR<=0.6 THEN XOGE*GFR/60 ELSE XOGE*GFR/3-0.19
; en plus
d/dt(ZCAE)=YCAI-YCA
init(ZCAE)=55
d/dt(ZMGE)=YMGI-YMG
init(ZMGE)=33
d/dt(ZSO4)=YSOI-YSO4
init(ZSO4)=11
d/dt(ZPO4)=YPOI-YPO4
init(ZPO4)=12.1; mM
;init(ZPO4)=21.78; mEq
d/dt(ZOGE)=YOGI-YORG
init(ZOGE)=66
d/dt(ZCLE)=YCLI-YCLU
init(ZCLE)=1144
XCAE=ZCAE/VEC
XMGE=ZMGE/VEC
XSO4=ZSO4/VEC
XPO4=ZPO4/VEC
XOGE=ZOGE/VEC
XCLE=ZCLE/VEC
XCLA=XCLE-DCLA
STBC=XCAE+XMGE-XSO4-1.8*XPO4-XOGE-XCLE+XNE+XKE-0.2214*XPP
;STBC=XCAE+XMGE-XSO4-XPO4-XOGE-XCLE+XNE+XKE-0.2214*XPP
YCLU=MAX(0,YNU+YKU-STPG+YNH4-YCO3+YCA+YMG-YSO4)
;*****************************************************************
; BLOCK 6 - KIDNEY 1
; Input = ADH, ALD, GFR, OSMP, PHA, THDF, XKE, XNE, YCO3, YGLU, YORG, YPO4, YMNU, YURU
; Output = STPG, QWU, YKU, YNH, YNH4,YNU
;*****************************************************************
; constant parameters
CPRX=0.2
YNH0=0.024
YTA0=0.0068
; F61
YNH4=YNH0*(-0.5*PHU1+4)
; F62
F62=YTA0*(-2.5*PHA1+19.5)
; F63
YTA1=IF PHU2<=4 THEN 0 ELSE IF (PHU2>4 AND PHU2<=5) THEN F62*(PHU2-4) ELSE F62
;en plus
d/dt(PHU2)=(PHU-PHU2)/TRSP
init(PHU2)=6
; F64
STPO=YPO4*(1+1/(1+10^(6.8-PHA)))
; F65 ;
GUESS PHU=6
ROOTI PHU=(YPO4*(10^-PHA+2*10^-6.8)/(10^-PHA+10^-6.8)-YPO4*(10^-PHU+2*10^-6.8)/(10^-PHU+10^-6.8))+(YORG*10^-4.3/(10^-PHA+10^-4.3)-YORG*10^-4.3/(10^-PHU+10^-4.3))-YTA
LIMIT PHU<=15
LIMIT PHU>=0
{A=STPG-YPO4
CORG=10^-4.3
CPO4=10^-6.8
B=A*(CPO4+CORG)-(CPO4*YPO4+CORG*YORG)
C=CPO4*CORG*(A-YPO4-YORG)
D=B^2-4*A*C
HU=(-B+D^0.5)/2/A
PHU=-log10(HU)}
{A=STPG-YPO4+CORG*YORG
CORG=1/(1+10^(PHA-4.3))
B=(STPG-YPO4+CORG*YORG)*(10^6.8+10^4.3)-YPO4*10^4.3-YORG*10^6.8
C=STPG-YPO4+CORG*YORG-YPO4-YORG
D=B^2-4*A*C
HU=(-B+D^0.5)/2/A
PHU=-log10(HU)}
; en plus
STPG=MAX(0,STPO+YORG-YTA)
YTA=YTA1+0.009+ALD*0.001
d/dt(PHA1)=(PHA-PHA1)/200
init(PHA1)=7.4
d/dt(PHU1)=(PHU-PHU1)/300
init(PHU1)=6
YNOD=MAX(0,YTA1+YNH4-YCO3) ;SRT- can't find YNOD in article or xls comparison file
YNU=MAX(YND*0.116-YNOD,0)
GFR1_6=THDF*GFR*CPRX
YND=XNE*GFR1_6*0.5*0.9-ALD*0.09
YKU=0.39*YKD
YKD=ALD*0.018*XKE+0.9*0.5*GFR1_6*XKE
OSMU=(YGLU+YURU+YMNU+1.86*(YKU+YNU))/QWU
QWD=(YGLU+YURU+YMNU+(YND+YKD)*1.86+0.32)/OSMP
QWU=QWD-QWD*0.9*ADH
YNH=0.5*GFR1_6*XNE
;*****************************************************************
; BLOCK 7 - CONTROLLER of RENAL FUNCTION
; Input = OSMP, PAS, PPCO, PVP, VEC, XKE, XNE,YNH
; Output = ADH, ALD, GFR, THDF
;*****************************************************************
; constant parameters
CKAL=0.5
CNAL=0.1
COAD=0.5
CPAD=1
CPAL=0.01
CPVL=0.1
GFR0=0.1
ACTH=1
VEC0=11
TADH=30; ou 15?
TALD=30
; F71
ADH=1.1/(1+exp(-0.5*(ADH0+4.605)))
; F72
ALD=10/(1+exp(-0.4394*(ALD0-5)))
; F73
THDF=IF PPCO<=28 THEN -5*(PPCO/28-1)+1 ELSE 1
; F74
GFR1=IF PAS<40 THEN 0 ELSE IF (PAS>=40 AND PAS <80) THEN 0.02*PAS-0.8 ELSE IF (PAS>=80 AND PAS<100) THEN -0.0005*(PAS-100)^2+1 ELSE 1
; en plus
GFR=GFR1*GFR0*VEC/VEC0
d/dt(ALD0)=(ALD1-ALD0)/TALD
init(ALD0)=0
d/dt(ADH0)=(z7-ADH0)/TADH
init(ADH0)=0
z7=(OSMP-287)*COAD-(PVP-4)*CPAD
;ALD1=(ACTH-1)*1+(XKE-4.5)*CKAL-(PVP-4)*CPVL-(YNH-1.4)*CNAL-(PAS-100)*CPAL
ALD1_av=(ACTH-1)*1+(XKE-4.5)*CKAL-(PVP-4)*CPVL-(YNH-1.4)*CNAL-(PAS-100)*CPAL
ALD1 = delay(ALD1_av, delay_time)
delay_time = 100
;figure12
;DISPLAY YGLI,XGLE, XKE,OSMP,QWU,YGLU,YKU
; figure11
;DISPLAY FCOI,FCOA,FO2A,PCOA,PO2A,VI,XCO3,UCOV,UCOA,PHA,PHU
;figure 10
;DISPLAY QVIN,ADH,ALD,STBC,PAS,VIF,OSMP,VIC,VEC,VP,QWU,QIN,PHU
