c f77 -o pa volca1.f implicit double precision (a-h,o-z) parameter(n=7) double precision y(n) double precision dydx(n), dym(n), dyt(n), yt(n) common /param/ vol, xica, rnaca, upt, dinf, zna, tauq common /p1/ cao, xnao, fz1,fz2, za, s12, xi, frt, xksat external derivs open(1,file='cp.dat') open(2,file='ci.dat') open(3,file='cj.dat') open(5,file='ica.dat') open(6,file='naca.dat') open(7,file='vol.dat') open(8,file='rel.dat') open(9,file='up.dat') open(11,file='f.dat') open(12,file='q.dat') open(13,file='d.dat') open(15,file='na.dat') open(16,file='imem.dat') open(17,file='cjp.dat') open(18,file='tq.dat') open(35,file='upsum.dat') open(36,file='exsum.dat') c------ runge stuff------------------------------------------ h=0.0001d0 ! time increment hh=h*0.5d0 ! used in runge algorithm h6=h/6.0d0 c------AP clamp --------------------------------------------- t0=0.10d0 tf=20.0d0 pcl=.40d0 ax=2.0d0/3.0d0 xr=ax/(ax+pcl) apd=xr*pcl volmin=-0.080d0 volmax=0.03d0 c-------sodium----------------------------------------------- znai=20.0d0 c znaf=13.6d0/(0.7d0+pcl) znaf=8.0d0 tauna=0.01d0 ! relaxation time of sodium c-------constants-------------------------------------------- cao=1.8d0 ! external ionic concentrations xnao=140.0d0 xf=96490.0d0 ! faraday's constant temp=35.0d0+273.0d0 ! temperature xr=8.315d0 ! gas constant R frt=xf/(xr*temp) 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-----------initial conditions---------------------------------------- y(1)=0.1d0 ! cp0 y(2)=0.10d0 ! ci0 y(3)=80.0d0 ! cj0 y(4)=.90d0 ! f starts at 1 y(5)=0.8d0 ! q y(6)=0.00d0 ! rate variable y(7)=y(3) ! jSR conc. x=t0 ! initialize time tt=t0 m=1 sx=0.0d0 ! initialize release sum do iter=1,100000000 tt=tt+h xm=m px1=xm*pcl ! begining of AP px2=px1+apd ! end of AP px3=px1+pcl if(tt.ge.px1.and.tt.le.px2) then ry=1.0d0-((tt-px1)/apd)**2 vol=(volmax-volmin)*sqrt(ry)+volmin else vol= volmin endif if(tt.ge.px1.and.tt.le.px3) then ! sub-membrane sodium zna=znaf+(-znaf+znai)*dexp(-(tt-px1)/tauna) else zna=znaf endif if(tt.gt.px3) then m=m+1 ! counter endif xh=x+hh if(tt.ge.tf) then ! end loop at final time goto 110 endif 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) pxdx1=px1+0.001d0 if(tt.gt.px1.and.tt.le.pxdx1) then upsum=0.0d0 exsum=0.0d0 casum=0.0d0 endif upsum=upsum+upt*h exsum=exsum+rnaca*h if(mod(iter,20).eq.0) then write(1,*) x , cp write(2,*) x , ci write(3,*) x , cj write(7,*) x , vol write(9,*) x , upt write(8,*) x , rel1 write(11,*) x , f write(12,*) x , q write(13,*) x , dinf write(14,*) x , uica write(15,*) x , zna write(16,*) x , rnaca-xica write(17,*) x, cjp write(18,*) x, tauq write(35,*) x , upsum write(36,*) x , exsum endif if(mod(iter,10).eq.0) then write(5,*) x , xica write(6,*) x , rnaca endif 100 end do 110 print*, 'finnish' print*, ' ' END c-----The derivatives are given here------------------------------ subroutine derivs(t,v,dvdt) implicit double precision (a-h,o-z) parameter(n=7) double precision t, v(n), dvdt(n) common /param/ vol, xica, rnaca, upt, dinf, zna, tauq common /p1/ cao, xnao, fz1,fz2, za, s12, xi, frt, xksat cp=v(1) ci=v(2) cj=v(3) f=v(4) q=v(5) rel=v(6) cjp=v(7) c snaca=0.50d0 snaca=4.0d0 sica=.50d0! scale Ica current for calcium entry part c------L-type channel current ------------------------------------- cai=ci/1000.0d0 cap=cp/1000.0d0 ! (mM) xnai=zna ! internal sodium conc. in mM volx=1000.0d0*vol ! d-gate vold=-5.0d0 dinf=1.0d0/(1.0d0+exp(-(volx+vold)/6.24d0)) ! d_inf ! f-gate volf=35.0d0 ! Luo-Rudy value finf=1.0d0/(1.0d0+exp((volx+volf)/8.6d0)) ! f_inf c tauf=0.030d0 ! slow inactivation in ms tauf=0.080d0 ! q-gate gam=1.0d0 cpstar=0.5d0 ! calcium inactivation threshold (micro M) c cpstar=1.50d0 qinf=1.0d0/(1.0d0+(cp/cpstar)**gam) c tauq=0.020d0 ! in seconds qb=8.0d0 tauq=qinf/qb za=vol*fz2 xis=(cap*exp(za)-0.341d0*cao)/(exp(za)-1.0d0) ! nernst current factor=vol*fz1 zica=factor*dinf*f*q*xis ! whole cell I_Ca current c------NaCa exchange---------------------------------------------- rt1=xi*vol se1=(rt1-vol)*frt ! exponents se2=rt1*frt se3=se1 s3=1.0d0+xksat*exp(se1) c s4=exp(se2)*xnai**3*cao-exp(se3)*xnao**3*cai 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------------------------------------------------------------------- fp=1.0d0/10.0d0 fi=1.0d0 fj=1.0d0/70.0d0 c------uptake------------------------------------------------------- vup=10.d0 xup=1.00d0 upt=vup*ci**2/(ci**2+xup**2) ! uptake c-----release------------------------------------------------------- taur=0.020d0 tauw=0.050d0 wx=1.50d0 gx=3.0d0 if(cjp.lt.50) then relj=0.0d0 else relj=wx*(cjp-50.0d0) endif cstar=110.0d0 av=10.0d0 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 c-----diffusion---------------------------------------------------- tau1=0.1d0 tr=(cp-ci)/tau1 ! cp-ci diffusion c-----I_Ca entry -------------------------------------------------- xica=sica*zica c-----leak--------------------------------------------------------- xleak=0.0d0 c-----final equations---------------------------------------------- gp=(rel - tr + rnaca - xica)/fp gi=(tr - upt + xleak )/fi gj=(-rel + upt - xleak)/fj gjp=(cj-cjp)/tauw ! jSR relaxation ff=(finf-f)/tauf ! f-gate fq=(qinf-q)/tauq ! q gate dvdt(1)=gp dvdt(2)=gi dvdt(3)=gj dvdt(4)=ff dvdt(5)=fq dvdt(6)=relx dvdt(7)=gjp return END