📄 step.java
字号:
/* This function performs the actual computations of currents, concentration changes
voltage changes, gate state changes involved in a single time step. Every 100 timesteps
the values of all currents are written to file. */
import java.lang.Math;
import java.io.File;
import java.io.*;
public class Step{
public static double knak;
public static double KmNa;
public static double KmK;
public static double knaca;
public static double KmNai;
public static double KmCa;
public static double ksat;
public static double n;
public static double Ko;
public static double Cao;
public static double Nao;
public static double Bufc;
public static double Kbufc;
public static double Bufsr;
public static double Kbufsr;
public static double taufca;
public static double taug;
public static double Vmaxup;
public static double Kup;
public static double pKNa;
public static double CAPACITANCE;
public static double RTONF;
public static double F;
public static double R;
public static double T;
public static double Gkr;
public static double Gks;
public static double GK1;
public static double Gto;
public static double GNa;
public static double GbNa;
public static double GCaL;
public static double GbCa;
public static double GpCa;
public static double KpCa;
public static double GpK;
public static double Vc;
public static double Vsr;
public static double period;
private double inverseVcF2;
private double inverseVcF;
private double Kupsquare;
private double BufcKbufc;
private double Kbufcsquare;
private double Kbufc2;
private double BufsrKbufsr;
private double Kbufsrsquare;
private double Kbufsr2;
private double exptaufca;
private double exptaug;
private double sm = 0;
private double sh = 0;
private double sj = 0;
private double sxr1 = 0;
private double sxr2 = 0;
private double sxs = 0;
private double ss = 0;
private double sr = 0;
private double sd = 0;
private double sf = 0;
private double sfca = 0;
private double sg = 0;
private double svolt= 0;
private double svolt2 = 0;
private double Cai = 0;
private double CaSR = 0;
private double Nai = 0;
private double Ki = 0;
private double sItot = 0;
private boolean flag = true;
public void step(Variable V,double HT,double tt,
int step,double Istim, FileOutputStream out)
{
//#define v(array_pointer,i,j) (*(V->array_pointer+i*V->NJ +j))
double IKr;
double IKs;
double IK1;
double Ito;
double INa;
double IbNa;
double ICaL;
double IbCa;
double INaCa;
double IpCa;
double IpK;
double INaK;
double Irel;
double Ileak;
double dNai;
double dKi;
double dCai;
double dCaSR;
double A;
double BufferFactorc;
double BufferFactorsr;
double SERCA;
double Caisquare;
double CaSRsquare;
double CaCurrent;
double CaSRCurrent;
double fcaold;
double gold;
double Ek;
double Ena;
double Eks;
double Eca;
double CaCSQN;
double bjsr;
double cjsr;
double CaBuf;
double bc;
double cc;
double Ak1;
double Bk1;
double rec_iK1;
double rec_ipK;
double rec_iNaK;
double AM;
double BM;
double AH_1;
double BH_1;
double AH_2;
double BH_2;
double AJ_1;
double BJ_1;
double AJ_2;
double BJ_2;
double M_INF;
double H_INF;
double J_INF;
double TAU_M;
double TAU_H;
double TAU_J;
double axr1;
double bxr1;
double axr2;
double bxr2;
double Xr1_INF;
double Xr2_INF;
double TAU_Xr1;
double TAU_Xr2;
double Axs;
double Bxs;
double Xs_INF;
double TAU_Xs;
double R_INF;
double TAU_R;
double S_INF;
double TAU_S;
double Ad;
double Bd;
double Cd;
double TAU_D;
double D_INF;
double TAU_F;
double F_INF;
double FCa_INF;
double G_INF;
// define all variables initialized
if (flag == true){
sm = V.M;
sh = V.H;
sj = V.J;
sxr1 = V.Xr1;
sxr2 = V.Xr2;
sxs = V.Xs;
ss = V.S;
sr = V.R;
sd = V.D;
sf = V.F;
sfca = V.FCa;
sg = V.G;
svolt = V.Volt;
svolt2 = V.Volt2;
Cai = V.Cai ;
CaSR = V.CaSR;
Nai = V.Nai;
Ki = V.Ki;
sItot = V.Itot;
flag = false;
}
inverseVcF2 = 1./(2*Vc*F);
inverseVcF = 1./(Vc*F);
Kupsquare = Kup*Kup;
BufcKbufc = Bufc*Kbufc;
Kbufcsquare = Kbufc*Kbufc;
Kbufc2 = 2*Kbufc;
BufsrKbufsr = Bufsr*Kbufsr;
Kbufsrsquare = Kbufsr*Kbufsr;
Kbufsr2 = 2*Kbufsr;
exptaufca = Math.exp(-HT/taufca);
exptaug = Math.exp(-HT/taug);
//Needed to compute currents
Ek=RTONF*(java.lang.Math.log((Ko/Ki)));
Ena=RTONF*(java.lang.Math.log((Nao/Nai)));
Eks=RTONF*(java.lang.Math.log((Ko+pKNa*Nao)/(Ki+pKNa*Nai)));
Eca=0.5*RTONF*(java.lang.Math.log((Cao/Cai)));
Ak1=0.1/(1.+Math.exp(0.06*(svolt-Ek-200)));
Bk1=(3.*Math.exp(0.0002*(svolt-Ek+100))+
Math.exp(0.1*(svolt-Ek-10)))/(1.+Math.exp(-0.5*(svolt-Ek)));
rec_iK1=Ak1/(Ak1+Bk1);
rec_iNaK=(1./(1.+0.1245*Math.exp(-0.1*svolt*F/(R*T))+0.0353*Math.exp(-svolt*F/(R*T))));
rec_ipK=1./(1.+Math.exp((25-svolt)/5.98));
//
//
//Compute currents
INa=GNa*sm*sm*sm*sh*sj*(svolt-Ena);
System.out.println(sm);
ICaL=GCaL*sd*sf*sfca*4*svolt*(F*F/(R*T))*
(Math.exp(2*svolt*F/(R*T))*Cai-0.341*Cao)/(Math.exp(2*svolt*F/(R*T))-1.);
Ito=Gto*sr*ss*(svolt-Ek);
IKr=Gkr*Math.sqrt(Ko/5.4)*sxr1*sxr2*(svolt-Ek);
IKs=Gks*sxs*sxs*(svolt-Eks);
IK1=GK1*rec_iK1*(svolt-Ek);
INaCa=knaca*(1./(KmNai*KmNai*KmNai+Nao*Nao*Nao))*(1./(KmCa+Cao))*
(1./(1+ksat*Math.exp((n-1)*svolt*F/(R*T))))*
(Math.exp(n*svolt*F/(R*T))*Nai*Nai*Nai*Cao-
Math.exp((n-1)*svolt*F/(R*T))*Nao*Nao*Nao*Cai*2.5);
INaK=knak*(Ko/(Ko+KmK))*(Nai/(Nai+KmNa))*rec_iNaK;
IpCa=GpCa*Cai/(KpCa+Cai);
IpK=GpK*rec_ipK*(svolt-Ek);
IbNa=GbNa*(svolt-Ena);
IbCa=GbCa*(svolt-Eca);
//
//
//Determine total current
sItot = IKr +
IKs +
IK1 +
Ito +
INa +
IbNa +
ICaL +
IbCa +
INaK +
INaCa +
IpCa +
IpK +
Istim;
//
//update concentrations
Caisquare=Cai*Cai;
CaSRsquare=CaSR*CaSR;
CaCurrent=-(ICaL+IbCa+IpCa-2*INaCa)*inverseVcF2*CAPACITANCE;
A=0.016464*CaSRsquare/(0.0625+CaSRsquare)+0.008232;
Irel=A*sd*sg;
Ileak=0.00008*(CaSR-Cai);
SERCA=Vmaxup/(1.+(Kupsquare/Caisquare));
CaSRCurrent=SERCA-Irel-Ileak;
CaCSQN=Bufsr*CaSR/(CaSR+Kbufsr);
dCaSR=HT*(Vc/Vsr)*CaSRCurrent;
bjsr=Bufsr-CaCSQN-dCaSR-CaSR+Kbufsr;
cjsr=Kbufsr*(CaCSQN+dCaSR+CaSR);
CaSR=(Math.sqrt(bjsr*bjsr+4*cjsr)-bjsr)/2;
CaBuf=Bufc*Cai/(Cai+Kbufc);
dCai=HT*(CaCurrent-CaSRCurrent);
bc=Bufc-CaBuf-dCai-Cai+Kbufc;
cc=Kbufc*(CaBuf+dCai+Cai);
Cai=(Math.sqrt(bc*bc+4*cc)-bc)/2;
dNai=-(INa+IbNa+3*INaK+3*INaCa)*inverseVcF*CAPACITANCE;
Nai+=HT*dNai;
dKi=-(Istim+IK1+Ito+IKr+IKs-2*INaK+IpK)*inverseVcF*CAPACITANCE;
Ki+=HT*dKi;
//
//
//write currents to file
if(step%100==0)
{
try {
String s ;
// s = "INa\t"+Double.toString(INa)+"\t"+
// "IKr\t"+Double.toString(IKr)+"\t"+
// "IKs\t"+Double.toString(IKs)+"\t"+
// "IK1\t"+Double.toString(IK1)+"\t"+
// "Ito\t"+Double.toString(Ito)+"\t"+
// "IbNa\t"+Double.toString(IbNa)+"\t"+
// "INaK\t"+Double.toString(INaK)+"\t"+
// "ICaL\t"+Double.toString(ICaL)+"\t"+
// "INaCa\t"+Double.toString(INaCa)+"\t"+
// "Irel\t"+Double.toString(Irel)+"\t"+
// "sItot\t"+Double.toString(sItot)+"\t"+
// "\n";
s = Double.toString(INa)+"\t"+
Double.toString(IKr)+"\t"+
Double.toString(IKs)+"\t"+
Double.toString(IK1)+"\t"+
Double.toString(Ito)+"\t"+
Double.toString(IbNa)+"\t"+
Double.toString(INaK)+"\t"+
Double.toString(ICaL)+"\t"+
Double.toString(INaCa)+"\t"+
Double.toString(Irel)+"\t"+
Double.toString(sItot)+"\t"+
"\n";
out.write(s.getBytes());
}catch(FileNotFoundException e){
}catch(IOException e){
}
}
//
//compute steady state values and time constants
AM=1./(1.+Math.exp((-60.-svolt)/5.));
BM=0.1/(1.+Math.exp((svolt+35.)/5.))+0.10/(1.+Math.exp((svolt-50.)/200.));
TAU_M=AM*BM;
M_INF=1./((1.+Math.exp((-56.86-svolt)/9.03))*(1.+Math.exp((-56.86-svolt)/9.03)));
if (svolt>=-40.){
AH_1=0.;
BH_1=(0.77/(0.13*(1.+Math.exp(-(svolt+10.66)/11.1))));
TAU_H= 1.0/(AH_1+BH_1);
} else{
AH_2=(0.057*Math.exp(-(svolt+80.)/6.8));
BH_2=(2.7*Math.exp(0.079*svolt)+(3.1e5)*Math.exp(0.3485*svolt));
TAU_H=1.0/(AH_2+BH_2);
}
H_INF=1./((1.+Math.exp((svolt+71.55)/7.43))*(1.+Math.exp((svolt+71.55)/7.43)));
if(svolt>=-40.)
{
AJ_1=0.;
BJ_1=(0.6*Math.exp((0.057)*svolt)/(1.+Math.exp(-0.1*(svolt+32.))));
TAU_J= 1.0/(AJ_1+BJ_1);
}
else
{
AJ_2=(((-2.5428e4)*Math.exp(0.2444*svolt)-(6.948e-6)*
Math.exp(-0.04391*svolt))*(svolt+37.78)/
(1.+Math.exp(0.311*(svolt+79.23))));
BJ_2=(0.02424*Math.exp(-0.01052*svolt)/(1.+Math.exp(-0.1378*(svolt+40.14))));
TAU_J= 1.0/(AJ_2+BJ_2);
}
J_INF=H_INF;
Xr1_INF=1./(1.+Math.exp((-26.-svolt)/7.));
axr1=450./(1.+Math.exp((-45.-svolt)/10.));
bxr1=6./(1.+Math.exp((svolt-(-30.))/11.5));
TAU_Xr1=axr1*bxr1;
Xr2_INF=1./(1.+Math.exp((svolt-(-88.))/24.));
axr2=3./(1.+Math.exp((-60.-svolt)/20.));
bxr2=1.12/(1.+Math.exp((svolt-60.)/20.));
TAU_Xr2=axr2*bxr2;
Xs_INF=1./(1.+Math.exp((-5.-svolt)/14.));
Axs=1100./(Math.sqrt(1.+Math.exp((-10.-svolt)/6)));
Bxs=1./(1.+Math.exp((svolt-60.)/20.));
TAU_Xs=Axs*Bxs;
// if( EPI)
R_INF=1./(1.+Math.exp((20-svolt)/6.));
S_INF=1./(1.+Math.exp((svolt+20)/5.));
TAU_R=9.5*Math.exp(-(svolt+40.)*(svolt+40.)/1800.)+0.8;
TAU_S=85.*Math.exp(-(svolt+45.)*(svolt+45.)/320.)+5./(1.+Math.exp((svolt-20.)/5.))+3.;
//#endif
//#ifdef ENDO
R_INF=1./(1.+Math.exp((20-svolt)/6.));
S_INF=1./(1.+Math.exp((svolt+28)/5.));
TAU_R=9.5*Math.exp(-(svolt+40.)*(svolt+40.)/1800.)+0.8;
TAU_S=1000.*Math.exp(-(svolt+67)*(svolt+67)/1000.)+8.;
//#endif
//#ifdef MCELL
R_INF=1./(1.+Math.exp((20-svolt)/6.));
S_INF=1./(1.+Math.exp((svolt+20)/5.));
TAU_R=9.5*Math.exp(-(svolt+40.)*(svolt+40.)/1800.)+0.8;
TAU_S=85.*Math.exp(-(svolt+45.)*(svolt+45.)/320.)+5./(1.+Math.exp((svolt-20.)/5.))+3.;
//#endif
D_INF=1./(1.+Math.exp((-5-svolt)/7.5));
Ad=1.4/(1.+Math.exp((-35-svolt)/13))+0.25;
Bd=1.4/(1.+Math.exp((svolt+5)/5));
Cd=1./(1.+Math.exp((50-svolt)/20));
TAU_D=Ad*Bd+Cd;
F_INF=1./(1.+Math.exp((svolt+20)/7));
TAU_F=1125*Math.exp(-(svolt+27)*(svolt+27)/240)+80+165/(1.+Math.exp((25-svolt)/10));
FCa_INF=(1./(1.+Math.pow((Cai/0.000325),8))+
0.1/(1.+Math.exp((Cai-0.0005)/0.0001))+
0.20/(1.+Math.exp((Cai-0.00075)/0.0008))+
0.23 )/1.46;
if(Cai<0.00035)
G_INF=1./(1.+Math.pow((Cai/0.00035),6));
else
G_INF=1./(1.+Math.pow((Cai/0.00035),16));
//Update gates
sm = M_INF-(M_INF-sm)*Math.exp(-HT/TAU_M);
sh = H_INF-(H_INF-sh)*Math.exp(-HT/TAU_H);
sj = J_INF-(J_INF-sj)*Math.exp(-HT/TAU_J);
sxr1 = Xr1_INF-(Xr1_INF-sxr1)*Math.exp(-HT/TAU_Xr1);
sxr2 = Xr2_INF-(Xr2_INF-sxr2)*Math.exp(-HT/TAU_Xr2);
sxs = Xs_INF-(Xs_INF-sxs)*Math.exp(-HT/TAU_Xs);
ss= S_INF-(S_INF-ss)*Math.exp(-HT/TAU_S);
sr= R_INF-(R_INF-sr)*Math.exp(-HT/TAU_R);
sd = D_INF-(D_INF-sd)*Math.exp(-HT/TAU_D);
sf =F_INF-(F_INF-sf)*Math.exp(-HT/TAU_F);
fcaold=sfca;
sfca =FCa_INF-(FCa_INF-sfca)*exptaufca;
if(sfca>fcaold && (svolt)>-60){
sfca=fcaold;
}
gold=sg;
sg =G_INF-(G_INF-sg)*exptaug;
if(sg>gold && (svolt)>-60){
sg=gold;
}
//update voltage
svolt= svolt + HT*(-sItot);
//System.out.println(svolt);
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -