c Compile with --- f77 -o pa cplca4z.f implicit double precision (a-h,o-z) parameter(n=15) parameter(npoints=1000000) parameter(kpoints=40) parameter(nbeat=70) double precision y(n), apdx(kpoints), dix(kpoints) double precision dydx(n), dym(n), dyt(n), yt(n) double precision cix(0:npoints), volxx(0:npoints), tx(0:npoints) double precision icax(0:npoints) double precision apd(0:nbeat), cimax(0:nbeat), di(0:nbeat) common /param/ xica, rnaca, upt, dinf, zna, tauq common /p1/ cao, xnao, fz1,fz2, za, s12, xi, frt, xksat common /p2/ ena, ek, eks, frtx, xkor, ckr, ystim, yina, yk common /p3/ yito, ykr, yks, yk1, ykp common /p4/ xbp, xbi, vito, vica, vicna, wica common /p5/ rxa external derivs c------output files------------------------------------------------- open(40,file='drest35.dat') open(41,file='cimax35.dat') open(42,file='brest35.dat') open(24,file='slope35.dat') c------ runge stuff------------------------------------------------ h=0.00005d0 ! time increment hh=h*0.5d0 ! used in runge algorithm h6=h/6.0d0 c------pacing stimulus ------------------------------------ volc=-0.070d0 ! threshold voltage to compute APD c-------internal sodium---------------------------------------------- zna=10.0d0 c-------constants for calcium system--------------------------------- cao=1.8d0 ! external ionic concentrations xnao=140.0d0 xnai=zna xf=96490.0d0 ! faraday's constant temp=35.0d0+273.0d0 ! temperature xr=8.315d0 ! gas constant R frt=xf/(xr*temp) frtx=frt/1000.0d0 ! convert to millivolts units c------parameters in I_Ca and I_NaCa--------------------------------- zs=2.0d0 ! charge pca=0.00054d0 fz1=pca*zs**2*xf**2/(xr*temp) fz2=zs*frt xkn=87.5d0 ! mM xkm=1.38d0 ! mM xi=0.35d0 xksat=0.1d0 s1=1.0d0/(xkn**3+xnao**3) s2=1.0d0/(xkm+cao) s12=s1*s2 c-------constants for voltage currents ------------------------------ xko=4.0d0 ! external pottassium milli molars xki=149.4d0 ! internal pottassium uc=1000.0d0*xr*temp/xf ! convert constant to millivolts ena=uc*Log(xnao/xnai) ! nernst potential for sodium in mV ek=uc*Log(xko/xki) eks=uc*Log((xko+0.01833d0*xnao)/(xki+0.01833*xnai)) xkmk1=13.0d0 xkor=xko/(xko+xkmk1) ckr=sqrt(xko/4.0d0) c-----------initial conditions---------------------------------------- ciox=0.4d0 cjox=110.0d0 y(1)=ciox ! cp0 ! calcium y(2)=ciox ! ci0 y(3)=cjox ! cj0 y(4)=0.90d0 ! f starts at 1 y(5)=0.7d0 ! q y(6)=0.00d0 ! rate variable y(7)=y(3) ! jSR conc. y(8)=-0.090d0 ! voltage y(9)= 0.99d0 ! sodium h-gate y(10)= 0.0d0 ! potassium xkr gate y(11)=0.0001d0 ! potassium xks gate y(12)=0.00001d0 ! potasssium xo gate y(13)=1.0d0 ! potassium yo gate y(14)=0.001d0 ! calcium d gate y(15)=1.0d0 ! sodum j gate c *************** dynamic protocol *************************************** pcli=1.20d0 c pcli=1.50d0 pclf=0.20d0 do kk=1, kpoints pcl=pcli+(pclf-pcli)*dfloat(kk-1)/dfloat(kpoints) c ********************************************************************** t0=0.0d0 x=t0 ! initialize time tt=t0 cix(0)=0.10d0 ! initialize apd(0)=0.100d0 cimax(0)=0.10d0 do m=1, nbeat ! beat iteration xm=dfloat(m) px1=xm*pcl ! begining of stimulus px2=px1+0.001d0 ! 1ms stimulus px3=px1+pcl do iter=1,100000000 ! within beat tt=tt+h if(tt.ge.px1.and.tt.le.px2) then ! stimulation current ystim=-50.0d0 ! stimulation current amplitude else ystim=0.0d0 endif xh=x+hh c----- Evolve using 4th order runge-kutta----------------- call derivs(x,y,dydx) do i=1,n yt(i)=y(i)+hh*dydx(i) end do call derivs(xh,yt,dyt) do i=1,n yt(i)=y(i)+hh*dyt(i) end do call derivs(xh,yt,dym) do i=1,n yt(i)=y(i)+h*dym(i) end do call derivs(xh,yt,dym) do i=1,n yt(i)=y(i)+h*dym(i) dym(i)=dyt(i)+dym(i) end do call derivs(x+h,yt,dyt) do i=1,n y(i)=y(i)+h6*(dydx(i)+dyt(i)+2.0d0*dym(i)) end do x=x+h ! increment the time cp=y(1) ci=y(2) cj=y(3) f=y(4) q=y(5) rel1=y(6) cjp=y(7) d=y(14) vol=y(8) cix(iter)=ci volxx(iter)=vol tx(iter)=tt icax(iter)=xica if(tt.ge.px3) then ! end beat loop at final time goto 50 endif end do 50 do k=3, iter-1 if(tx(k).gt.(0.02)) then if(volxx(k+1).lt.volc.and.volxx(k).gt.volc) then kapd=k endif endif if(cix(k).gt.cix(k-1).and.cix(k).gt.cix(k+1)) then kmax=k endif end do vx1=volxx(kapd) vx2=volxx(kapd+1) tc=((volc-vx1)*h)/(vx2-vx1) ! interpolate apd(m)=tx(kapd) + tc - px1 cimax(m)=cix(kmax) di(m)=pcl-apd(m) if(volxx(iter-1).gt.volc) then apd(m)=0.0d0 di(m)=0.0d0 endif if(m.gt.(nbeat-2)) then ! last two beats write(40,*) di(m-1), apd(m) write(41,*) pcl, cimax(m) write(42,*) pcl, apd(m) endif end do ! goto next beat apdx(kk)=apd(nbeat) ! apd and di at last beat dix(kk)=di(nbeat-1) end do ! compute dynamic restitution slope here do kk=2, kpoints c if(dix(kk).gt.(0.015)) then pcl=pcli+(pclf-pcli)*dfloat(kk-1)/dfloat(kpoints) slope=(apdx(kk)-apdx(kk-1))/(dix(kk)-dix(kk-1)) write(24,*) pcl, slope c endif end do END c-----The derivatives are given here------------------------------ subroutine derivs(t,v,dvdt) implicit double precision (a-h,o-z) parameter(n=15) double precision t, v(n), dvdt(n) common /param/ xica, rnaca, upt, dinf, zna, tauq common /p1/ cao, xnao, fz1,fz2, za, s12, xi, frt, xksat common /p2/ ena, ek, eks, frtx, xkor, ckr, ystim, yina, yk common /p3/ yito, ykr, yks, yk1, ykp common /p4/ xbp, xbi, vito, vica, vicna, wica common /p5/ rxa cp=v(1) ci=v(2) cj=v(3) f=v(4) q=v(5) rel=v(6) cjp=v(7) vol=v(8) h=v(9) ! sodium h-gate xkr=v(10) xks=v(11) xo=v(12) yo=v(13) d=v(14) ! calcium d gate xj=v(15) snaca=4.0d0 sica=3.50d0 gxr=0.30d0 ! scale g of release separately wca=0.02d0 ! controls stength of calcium currents wk=1.00d0 ! adjustable constant to control AP shape; this controls ! the strength of the potassium currents c******************************************************************** volx=1000.0d0*vol ! convert from volts to millivolts c**************** L-type channel current ************************************ cap=cp/1000.0d0 ! (mM) xnai=zna ! internal sodium conc. in mM c--------- d-gate ------------------------------------------------- vold=-5.0d0 dinf=1.0d0/(1.0d0+exp(-(volx+vold)/6.24d0)) ! d_inf taud=0.005d0 ! fast activation of the d-gate c---------f-gate--------------------------------------------------- volf=35.0d0 ! Luo-Rudy value finf=1.0d0/(1.0d0+dexp((volx+volf)/8.6d0)) ! f_inf tauf=0.035d0 c---------q-gate--------------------------------------------------- gam=2.0d0 cpstar=0.50d0 qinf=1.0d0/(1.0d0+(cp/cpstar)**gam) tauq=0.02d0 c------------------------------------------------------------------- za=vol*fz2 factor=vol*fz1 rxa=factor*(cap*dexp(za)-0.341d0*cao)/(dexp(za)-1.0d0) ! nernst current zica=d*f*q*rxa ! whole cell I_Ca current c zica=d*f*rxa c------NaCa exchange------------------------------------------------ rt1=xi*vol se1=(rt1-vol)*frt ! exponents se2=rt1*frt se3=se1 s3=1.0d0+xksat*exp(se1) s4=exp(se2)*xnai**3*cao-exp(se3)*xnao**3*cap xnaca=2000.0d0 ! can be adjusted rnaca=snaca*xnaca*s12*s4/s3 ! NaCa exchange current c------uptake------------------------------------------------------- vup=250.0d0 ! uptake strength xup=0.50d0 ! uptake treshold upt=vup*ci**2/(ci**2+xup**2) ! uptake current c------The SR release function------------------------------------- taur=0.020d0 ! spark lifetime tauw=0.050d0 ! relaxation time between jSR and SR wx=1.500d0 ! release function constant gx=gxr*250.0d0 ! release function constant if(cjp.lt.50) then relj=0.0d0 else relj=wx*(cjp-50.0d0) endif cstar=110.0d0 c av=30.0d0 av=2.0d0 ! stable calcium bv=wx*(cstar-50.0d0)-av*cstar if(cjp.ge.cstar) then relj=av*cjp+bv endif rr1=gx*(-zica)*relj relx=rr1-rel/taur ! ode for release flux if(cj.lt.10.0d0) then rel=0.0d0 endif c------Diffusion from submembrane to myoplasm----------------------- tau1=0.002d0 tr=(cp-ci)/tau1 ! cp-ci diffusion c xxxxxxxxxxxxxxx I_Ca entry current xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xica=sica*zica wica=xica/q c xxxxxxxxxxxxxxxxx buffers xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx bcal=24.0d0 ! calmodulin in the cytosol xkcal=7.0d0 srmax=47.0d0 srkd=0.60d0 timax=70.0d0 tkd=0.60d0 bpx=bcal*xkcal/(xkcal+cp)**2 ! calmodulin spx=srmax*srkd/(srkd+cp)**2 ! add SR buffering tpx=timax*tkd/(tkd+cp)**2 ! instantaneous Troponin bp=1.0d0/(1.0d0+bpx+spx+tpx) bix=bcal*xkcal/(xkcal+ci)**2 six=srmax*srkd/(srkd+ci)**2 ! add SR buffering tix=timax*tkd/(tkd+ci)**2 ! instantaneous Troponin bi=1.0d0/(1.0d0+bix+six+tix) vis=10.0d0 ! volume ratio vi/vs gp=bp*( vis*(rel - tr + rnaca - xica)) gi=bi*(tr - upt) gj=-rel + upt gjp=(cj-cjp)/tauw ! jSR relaxation ff=(finf-f)/tauf ! f-gate fq=(qinf-q)/tauq ! q-gate gd=(dinf-d)/taud ! d-gate c xxxxxxxxxxxxxSodium currentxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ! h gate alph=0.135d0*exp((volx+80.0d0)/(-6.8d0)) beth=7.5d0/(1.0d0+exp(-0.1d0*(volx+11.0d0))) hinf=alph/(alph+beth) tauh1=1.0d0/(alph+beth) tauh=tauh1/1000.0d0 ! convert to seconds gh=(hinf-h)/tauh ! ode for sodium h-gate c---------------------------------------------------------------- ! m-gate m=m_infinity sr1=-0.1d0*(volx+47.13d0) alpm=0.32d0*(volx+47.13d0)/(1.0d0-exp(sr1)) betm=0.08d0*exp(-volx/11.0d0) zminf=alpm/(alpm+betm) c--------------------------------------------------------------- ! j gate rs1=(volx+100.0d0)/(-23.0d0) rs2=0.15d0*(volx+79.0d0) alpj=0.175d0*exp(rs1)/(1.0d0+exp(rs2)) rs3=-0.1d0*(volx+79.0d0) betaj=0.3d0/(1.0d0+exp(rs3)) xjinf=alpj/(alpj+betaj) tauj1=1.0d0/(alpj+betaj) tauj=tauj1/1000.0d0 gxj=(xjinf-xj)/tauj c--------------------------------------------------------------- c gna=12.80d0 c gna=10.0d0 gna=12.0d0 yina=gna*(zminf**3)*h*xj*(volx-ena) ! sodium current c xxxxxxxxxxxx Pottasium currents xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ! Time independent potassium current I_K1 yz1=0.06d0*(volx-ek) yk1inf=1.0d0/(2.0d0+exp(yz1)) gk1=2.8d0 c gk1=3.0d0 yk1=gk1*yk1inf*xkor*(volx-ek) ! Time dependent potassium current I_Kr xkrinf=1.0d0/(1.0d0+exp(-2.182d0-0.1819*volx)) zr5=exp(-5.495+0.1691*volx)+exp(-7.677-0.0128*volx) taukr1=43.0d0+1.0d0/zr5 taukr=taukr1/1000.0d0 ! convert to seconds gxkr=(xkrinf-xkr)/taukr ! ode for xkr gate xrv=1.0d0/(1.0d0+2.5d0*exp(0.10d0*(volx+28.0d0))) gkr=0.0136d0 c gkr=0.06d0 ykr=gkr*xrv*xkr*ckr*(volx-ek) ! Plateau potassium current gkp=0.002216d0 c gkp=0.001d0/4.0d0 xkp=1.0d0/(1.0d0+exp((7.488d0-volx)/5.98d0)) ykp=gkp*xkp*(volx-ek) c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ! I_Ks gks=0.0245d0 c gks=0.01d0 ar1=(volx-16.0d0)/(-13.6d0) xksinf=1.0d0/(1.0d0+exp(ar1)) as2=0.0000719d0*(volx-10.0)/(1.0d0-exp(-0.148d0*(volx-10.0))) as3=0.000131d0*(volx-10.0)/(exp(0.0687*(volx-10.0))-1.0d0) as1=as2+as3 tauks1=1.0d0/as1 tauks=tauks1/1000.0d0 ! convert to seconds gxks=(xksinf-xks)/tauks yks=gks*xks**2*(volx-eks) c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ! I_to c gito=0.23815d0 c gito=0.07d0 gito=0.10d0 axto=0.04516d0*exp(0.03577*volx) bxto=0.0989*exp(-0.06237d0*volx) ar1=(volx+33.5)/(-5.) ayto1=0.005415d0*exp(ar1) ayto2=1.0d0+0.051335*exp(ar1) ayto=ayto1/ayto2 byto1=0.005415d0*exp(-ar1) byto2=1.0d0+0.051335*exp(-ar1) byto=byto1/byto2 xoinf=axto/(axto+bxto) txo1=1.0d0/(axto+bxto) txo=txo1/1000.0d0 ! convert to seconds yoinf=ayto/(ayto+byto) tyo1=1.0d0/(ayto+byto) tyo=tyo1/1000.0d0 ! convert to seconds gxto=(xoinf-xo)/txo gyto=(yoinf-yo)/tyo yito=gito*xo*yo*(volx-ek) ykw=yk1+ykr+ykp+yks+yito ! total potassium current c xxxxxxxxxxxxxx Equations for voltage xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx yca=wca*(2.0d0*xica+rnaca) ! total calcium current yk=wk*ykw ! total potassium current gvol=-(yina+yk+yca+ystim) vito=wk*yito vica=wca*xica vicna=wca*(xica+rnaca) dvdt(1)=gp dvdt(2)=gi dvdt(3)=gj dvdt(4)=ff dvdt(5)=fq dvdt(6)=relx dvdt(7)=gjp dvdt(8)=gvol dvdt(9)=gh dvdt(10)=gxkr dvdt(11)=gxks dvdt(12)=gxto dvdt(13)=gyto dvdt(14)=gd dvdt(15)=gxj return END