- Author:
- Randall Britten <r.britten@auckland.ac.nz>
- Date:
- 2014-03-10 12:37:50+13:00
- Desc:
- Snapshot of partial progress towards a CellML representation, based on the BM file.
At this stage, this aimed just at the Respiration component, and without physical units, or any work to validate that the simulation results are correct.
A more detailed version history to this point is available at https://github.com/codecurve/Ikeda-et-al-1979
- Permanent Source URI:
- https://models.cellml.org/workspace/193/rawfile/e837300b086b879f3016fa777f83b9199afd1ee8/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
