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