Biophysically detailed model of the mouse sino-atrial node cell (Kharche et al. 2011)

 Download zip file 
Help downloading and running models
Accession:141274
This model is developed to study the role of various electrophysiological mechanisms in generating cardiac pacemaking action potentials (APs).The model incorporates membrane ionic currents and intracellular mechanisms contributing to spontaneous mouse SAN APs. The model was validated by testing the functional roles of individual membrane currents in one and multiple parameter analyses.The roles of intracellular Ca2+-handling mechanisms on cardiac pacemaking were also investigated in the model.
Reference:
1 . Kharche S, Yu J, Lei M, Zhang H (2011) A mathematical model of action potentials of mouse sinoatrial node cells with molecular bases. Am J Physiol Heart Circ Physiol 301:H945-63 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Cardiac atrial cell;
Channel(s): I N; I K; I Sodium; I Calcium; I Potassium; I_HERG; Na/Ca exchanger; KCNQ1; Na/K pump;
Gap Junctions:
Receptor(s):
Gene(s): Nav1.1 SCN1A; Nav1.2 SCN2A; HERG KCNH2; Cav3.1 CACNA1G; Cav1.3 CACNA1D; Cav1.2 CACNA1C;
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Action Potential Initiation; Oscillations; Action Potentials; Calcium dynamics; Cardiac pacemaking;
Implementer(s): Kharche, Sanjay ;
Search NeuronDB for information about:  I N; I K; I Sodium; I Calcium; I Potassium; I_HERG; Na/Ca exchanger; KCNQ1; Na/K pump;
/
Kharche_SAN
readme.html
Kharche_SAN.c
screenshot.jpg
screenshot1.jpg
                            
/*
July 2011.
A mathematical model of action potentials of mouse sinoatrial node cells with molecular bases
Sanjay Kharche, Jian Yu, Ming Lei, and Henggui Zhang.
AJP - Heart September 2011 vol. 301 no. 3 H945-H963

LINK TO PAPER:

http://ajpheart.physiology.org/content/301/3/H945.long
*/

#include <stdio.h>
#include <math.h>
#include<stdlib.h>

#define R 8.314472
#define T 310.5   
#define F 96.4845

#define number_of_apds 100

/***************************************************************************************/

#define ddt 0.00025
#define capacitance 0.025
#define vcell 3.0
#define l_cell 66.3767257
#define r_cell 3.792956
#define vrel 0.0036
#define vsub 0.03328117
#define vup  0.0348
#define vi 1.34671883
#define Mgi 2.5
#define nao 140.0
#define cao 1.8
#define ko 5.4
#define gst 0.00006
#define eist 17.0
#define gbna 0.0001215
#define gbca 0.000015
#define gbk  0.0000025
#define	gk1 0.229*0.0039228*0.9
#define gks 0.000299
#define ecal 47.0
#define kmfca 0.00035
#define alpha_fca 0.021
#define all_ica_multiplier 1.0
#define ecat 45.0
#define enattxr  41.5761
#define multiplier2 1.0
#define gsus 0.00039060
#define inakmax_multiplier 1.85
#define inakmax inakmax_multiplier*0.077
#define kmnap 14.0
#define kmkp 1.4
#define K1ni 395.3
#define K1no 1628
#define K2ni 2.289
#define K2no 561.4
#define K3ni 26.44
#define K3no 4.663
#define Kci 0.0207
#define Kco 3.663
#define Kcni 26.44
#define Qci 0.1369
#define Qco 0.0
#define Qn 0.4315
#define tdifca 0.04
#define Prel 2.5
#define Krel 0.0015
#define nrel 2.0
#define Kup 0.0006
#define nup 1.0
#define Ttr 40.0
#define ConcTC 0.031
#define ConcTMC 0.062
#define kfTC 88.8
#define kfTMC 237.7
#define kbTC 0.446
#define kbTMC 0.00751
#define kfTMM 2.277
#define kbTMM 0.751
#define ConcCM 0.045
#define kfCM 237.7
#define kbCM 0.542
#define ConcCQ 10.0
#define kfCQ 0.534
#define kbCQ 0.445
#define koca 10.0
#define kom 0.06
#define kica 0.5
#define kim 0.005
#define eca50sr 0.45
#define maxsr 15.0
#define minsr 1.0
#define hsrr 2.5
#define pumphill 2.0

double v, vnew,sstime;
int  current_trace_counter = 0, output_counter = 0;
double vclamp;
double total_current;
double cai; 
double casub;
double ena,eca,ek,eks;
double ist;
double qa, qi, tauqa,tauqi,alphaqa,betaqa,alphaqi,betaqi;
double dst;
double fst;
double ibca,ibna,ibk,ib;
double ik1 ;
double xk1inf;
double icat;
double dt;
double ft;
double dt_inf, tau_dt, ft_inf, tau_ft;
double ikr_act, ikr_act_inf,tau_ikr_act;
double ikr;
double ikr_inact_inf,tau_ikr_inact,tau_ikr_inact2,ikr_inact,ikr_inact2;
double iks_act, iks_act_inf,tau_iks_act;
double iks;
double alpha_dl, beta_dl, tau_dl, alpha_fl, beta_fl, tau_fl;
double fl12, dl12, dl12_inf, fl12_inf,ical12;
double fl13, dl13, dl13_inf, fl13_inf,ical13;
double fca; 
double fca_inf;
double taufca;
double ina_ttxr, ina_ttxs;
double m3_inf_ttxr, h_inf_ttxr;
double m3_inf_ttxs, h_inf_ttxs;
double m_inf_ttxr,m_inf_ttxs, hs,hsr;
double m_ttxr,h_ttxr,j_ttxr;
double m_ttxs,h_ttxs,j_ttxs;
double tau_m,tau_h,tau_j,tau_mr,tau_hr,tau_jr;
double delta_ttxr_m,delta_ttxr_h,delta_ttxr_j;
double fna;
double ih, ihk, ihna, ih_1, ih_2, ih_4, y_1_2, y_4, tau_y_1_2, tau_y_4, y_inf, y_1_2_inf, y_4_inf,ih_1_k, ih_2_k, ih_4_k, ih_1_na, ih_2_na, ih_4_na;
double inak;
double inaca;
double di,doo,k43,k12,k14,k41,k34,k21,k23,k32,x1,x2,x3,x4;
double r_inf, r, tau_r, isus;
double q,q_inf,tau_q;
double ito;
double ca_flux;
double Jcadif;
double Jrel, Jup, Jtr,carel,caup;
double variation_multiply = 1.0;
double dFtc,Ftc;
double dFtmc,Ftmc,Ftmm;
double dFtmm,Ftmm;
double dFcms,Fcms;
double dFcmi,Fcmi;
double dFcq,Fcq;
int param_counter = 0, number = 0;
double min_potential[number_of_apds];
double tmin_potential[number_of_apds];
double max_potential[number_of_apds];
double dvdtmax[number_of_apds];
double vdvdtmax[number_of_apds];
double apd_start[number_of_apds], apd_end[number_of_apds],cai_peak[number_of_apds],cai_min[number_of_apds];
double ddr[number_of_apds];
double top[number_of_apds];
double top_slope[number_of_apds];
double apd50[number_of_apds];
double apd90[number_of_apds];
double cycle_length[number_of_apds];
double casub_peak[number_of_apds],casub_min[number_of_apds];
double dummy = 1.0;
double base_cycle = 0.0;
double dvdt, dvdtold;
double nai_tot, ki_tot;
double nai_integral, ki_integral;
double ki_integral1, ki_integral2;
double nai_min, nai_max, ki_min, ki_max;
double nai=8.0;  
double ki=140;
double FRT, RTF;
double resting, open, inactivated, resting_inactivated;
double kcasr,kosrca,kisrca;
double pumpvmax = 0.04;
double pumpkmf = 0.000246;
double pumpkmr = 3.29;
double pumpfactor = 1.0;
double Pup;
double Pup2;
double ks;
double gna_ttxs;
double gna_ttxr;
double gcat;
double gkr;
double gto;
double gcal12;
double gcal13;
double gh;
double vhalf_gh;
double kNaCa;

int main()
{  
FILE *outputcurrents,*parameters, *random;

char *str;
int start_output = 0;
double total_time = 0.0;
int temp_counter = 0;

FRT = F/(R*T);
param_counter = 0;
current_trace_counter = 0;
outputcurrents = fopen("mousesan.dat","w+");

/*
Initial conditions
*/

v = -64.5216286940;
dst = 0.6246780312;
fst = 0.4537033169;
dt = 0.0016256324;
ft = 0.4264459666;
ikr_act = 0.4043600437;
ikr_inact = 0.9250035423;
ikr_inact2 = 0.1875749806;
iks_act = 0.0127086259;
fl12 = 0.9968141226;
dl12 = 0.0000045583;
fl13 = 0.9809298233;
dl13 = 0.0002036671;
r = 0.0046263658;
m_ttxr = 0.4014088304;
h_ttxr = 0.2724817537;
j_ttxr = 0.0249208708;
m_ttxs = 0.1079085266;
h_ttxs = 0.4500098710;
j_ttxs = 0.0268486392;
y_1_2 = 0.0279984462;
y_4 = 0.0137659036;
carel = 0.1187281829;
caup = 1.5768287365;
casub = 0.0000560497;
Ftc = 0.0063427103;
Ftmc = 0.1296677919;
Ftmm = 0.7688656371;
Fcms = 0.0242054739;
Fcmi = 0.0138533048;
Fcq = 0.1203184861;
cai = 0.0000319121;
q = 0.6107148187;
fca = 0.7649576191;
nai = 8.1179761505;
ki = 139.8854603066;
resting = 0.7720290515;
open = 0.0000000760;
inactivated = 0.0000000213;
resting_inactivated = 0.2162168926;

gna_ttxs = 0.1*5.925e-05;
gna_ttxr = 0.1*5.925e-05;
gcal12 = 0.0010*4.0*1.5;
gcal13 = 0.0030*4.0*1.5;
gcat = 0.75*0.01862; 
gh = 0.0057; 
vhalf_gh = 106.8;
gkr = 0.8*0.002955;
gto = 0.00492;
 kNaCa = 5.5;
Pup = 0.04;
ks = 1300000.0;
Pup2 = 0.04;
pumpkmf = 0.00008;
pumpkmr = 4.5;

for(param_counter=0;param_counter<number_of_apds;param_counter++){
     min_potential[param_counter] = 100000.0;
     max_potential[param_counter] = -10000.0;
     dvdtmax[param_counter]       = -10000.0;
     ddr[param_counter]           = -10000.0;
     top[param_counter]           =  10000.0;
     top_slope[param_counter]     = -10000.0;
     apd50[param_counter]         = -10000.0;
     apd90[param_counter]         = -10000.0;
     cycle_length[param_counter]  = -10000.0;
    }

ena = (R*T/F)*log(nao/nai);
ek  = (R*T/F)*log(ko/ki);
eks = ((R*T)/F)*log((ko+0.12*nao)/(ki+0.12*nai));
eca = (R*T/(2*F))*log(cao/casub);

param_counter = 0;
dvdt = 1000.0;
dvdtold = 500.0;
start_output = 0;
total_time = 5000.0;

for(sstime=0.0;sstime<=total_time;sstime = sstime + ddt){

/*Ist********************************************************************/

qa = 1.0/(1.0 + exp(-(v+67.0)/5.0));
alphaqa = 1.0/(0.15*exp(-(v)/11.0)+0.2*exp(-(v)/700.0));
betaqa  =  1.0/(16.0*exp((v)/8.0)+15.0*exp((v)/50.0));
tauqa = 1.0/(alphaqa + betaqa);
alphaqi = 0.15*1.0/(3100.0*exp((v+10.0)/13.0)+700.3*exp((v+10.0)/70.0));
betaqi =  0.15*1.0/(95.7*exp(-(v+10.0)/10.0) + 50.0*exp(-(v+10.0)/700.0)) + 0.000229/(1+exp(-(v+10.0)/5.0));
qi = alphaqi/(alphaqi + betaqi);
tauqi = 1.0/(alphaqi + betaqi);
dst = dst + ddt*((qa-dst)/tauqa);
fst = fst + ddt*((qi-fst)/tauqi);
ist = gst*dst*fst*(v - eist);

/* Ib ************************************************************************/

ibna = gbna*(v - ena);
ibca = gbca*(v - eca);
ibk  =  gbk*(v - ek);
ib = (ibna + ibca + ibk);
  
/*IK1**********************************************************************/

xk1inf = 1.0/(1.0 + exp(0.070727*(v - ek)));
ik1 = gk1*xk1inf*(ko/(ko + 0.228880))*(v - ek);

/**ICaT Cav3.1**************************************************************/

tau_dt = 1.0/(1.068*exp((v + 26.3)/30.0) + 1.068*exp(-(v + 26.3)/30.0));
dt_inf = 1.0/(1.0+exp(-(v + 26.0)/6.0));
dt = dt + ddt*((dt_inf - dt)/tau_dt);
tau_ft = 1.0/(0.0153*exp(-(v+61.7)/83.3)+0.015*exp((v+61.7)/15.38));
ft_inf = 1.0/(1.0+exp((v + 61.7)/5.6));
ft = ft + ddt*((ft_inf - ft)/tau_ft);
icat = gcat*ft*dt*(v - ecat);          

/*Ikr********************************************************************/

ikr_act_inf = 1.0/(1.0 + exp(-(v+21.173694)/9.757086));
tau_ikr_act = 0.699821/(0.003596*exp((v)/15.339290) + 0.000177*exp(-(v)/25.868423));
ikr_act = ikr_act + ddt*(ikr_act_inf-ikr_act)/tau_ikr_act;     
ikr_inact_inf = 1.0/(1.0 + exp((v+20.758474-4.0)/(19.0)));
tau_ikr_inact = 0.2+0.9*1.0/(0.1*exp(v/54.645)+0.656*exp(v/106.157));
ikr_inact = ikr_inact + ddt*(ikr_inact_inf - ikr_inact)/tau_ikr_inact;
ikr = gkr*ikr_act*ikr_inact*(v - ek);

/**IKs********************************************************************/

iks_act_inf = 1.0/(1.0 + exp(-(v-20.876040)/11.852723));
tau_iks_act =  1000.0/(13.097938/(1.0 + exp(-(v-48.910584)/10.630272)) + exp(-(v)/35.316539));
iks_act = iks_act + ddt*(iks_act_inf - iks_act)/tau_iks_act;
iks = gks*iks_act*iks_act*(v - eks);

/*ICaL*******************************************************************/

if(fabs(v)<=0.001){
alpha_dl  = -28.39*(v+35.0)/(exp(-(v+35.0)/2.5)-1.0)+408.173;
}
else
if(fabs(v+35.0)<=0.001){
alpha_dl  = 70.975-84.9*v/(exp(-0.208*v)-1.0);
}
else
if(fabs(v)>0.001&&fabs(v+35.0)>0.001){
alpha_dl  = -28.39*(v+35.0)/(exp(-(v+35.0)/2.5)-1.0)-84.9*v/(exp(-0.208*v)-1.0);
}
if(fabs(v-5.0)<=0.001)
 beta_dl   = 28.575;
else
if(fabs(v-5.0)>0.001)
beta_dl   = 11.43*(v-5.0)/(exp(0.4*(v-5.0))-1.0);
tau_dl  = 2000.0/(alpha_dl +beta_dl);
dl13_inf = 1.0/(1+exp(-(v+13.5)/6.0));
fl13_inf = 1.0/(1+exp((v+35.0)/7.3));
tau_fl = (7.4 + 45.77*exp(-0.5*(v+28.1)*(v+28.1)/(11*11)));
dl13 = dl13 + ddt*(dl13_inf - dl13)/tau_dl;
fl13 = fl13 + ddt*(fl13_inf - fl13)/tau_fl;
dl12_inf = 1.0/(1+exp(-(v+3.0)/5.0));
fl12_inf = 1.0/(1+exp((v+36.0)/4.6));
dl12 = dl12 + ddt*(dl12_inf - dl12)/tau_dl;
fl12 = fl12 + ddt*(fl12_inf - fl12)/tau_fl;
fca_inf = kmfca/(kmfca+casub);
taufca = fca_inf/alpha_fca;
fca = fca + ddt*(fca_inf - fca)/taufca;
ical12 = gcal12*fl12*dl12*fca*(v-ecal);
ical13 = gcal13*fl13*dl13*fca*(v-ecal);
   
/**INa**********************************************************************/

fna = (9.52e-02*exp(-6.3e-2*(v+34.4))/(1+1.66*exp(-0.225*(v+63.7))))+8.69e-2; 
m3_inf_ttxr = 1.0/(1.0 + exp(-(v+45.213705)/7.219547));
h_inf_ttxr = 1.0/(1.0 + exp(-(v+62.578120 )/(-6.084036)));
m3_inf_ttxs = 1.0/(1.0 + exp(-(v+36.097331-5.0)/5.0));
h_inf_ttxs = 1.0/(1.0 + exp((v+56.0)/3.0));
m_inf_ttxr = pow(m3_inf_ttxr,0.333);
m_inf_ttxs = pow(m3_inf_ttxs,0.333);
tau_m = 1000.0*((0.6247e-03/(0.832*exp(-0.335*(v+56.7))+0.627*exp(0.082*(v+65.01))))+0.0000492);
tau_h = 1000.0*(((3.717e-06*exp(-0.2815*(v+17.11)))/(1+0.003732*exp(-0.3426*(v + 37.76))))+0.0005977);
tau_j = 1000.0*(((0.00000003186*exp(-0.6219*(v+18.8)))/(1+0.00007189*exp(-0.6683*(v+34.07))))+0.003556);
m_ttxs = m_ttxs + ddt*(m_inf_ttxs - m_ttxs)/tau_m;
h_ttxs = h_ttxs + ddt*(h_inf_ttxs - h_ttxs)/tau_h;
j_ttxs = j_ttxs + ddt*(h_inf_ttxs - j_ttxs)/tau_j;
hs = (1.0-fna)*h_ttxs+fna*j_ttxs;
tau_mr = 1000.0*((0.6247e-03/(0.832*exp(-0.335*(v+56.7))+0.627*exp(0.082*(v+65.01))))+0.0000492);
tau_hr = 1000.0*(((3.717e-06*exp(-0.2815*(v+17.11)))/(1+0.003732*exp(-0.3426*(v + 37.76))))+0.0005977);
tau_jr = 1000.0*(((0.00000003186*exp(-0.6219*(v+18.8)))/(1+0.00007189*exp(-0.6683*(v+34.07))))+0.003556);
m_ttxr = m_ttxr + ddt*(m_inf_ttxr - m_ttxr)/tau_mr;
h_ttxr = h_ttxr + ddt*(h_inf_ttxr - h_ttxr)/tau_hr;
j_ttxr = j_ttxr + ddt*(h_inf_ttxr - j_ttxr)/tau_jr;
hsr = (1.0-fna)*h_ttxr+fna*j_ttxr;
if(fabs(v)>0.005)
ina_ttxs= gna_ttxs*m_ttxs*m_ttxs*m_ttxs*hs*nao*(F*F/(R*T))*((exp((v-ena)*F/(R*T))-1.0)/(exp(v*F/(R*T))-1.0))*v;
else
ina_ttxs= gna_ttxs*m_ttxs*m_ttxs*m_ttxs*hs*nao*F*((exp((v-ena)*F/(R*T))-1.0));
if(fabs(v)>0.005)
ina_ttxr = gna_ttxr*m_ttxr*m_ttxr*m_ttxr*hsr*nao*(F*F/(R*T))*((exp((v-enattxr)*F/(R*T))-1.0)/(exp(v*F/(R*T))-1.0))*v;
else
ina_ttxr = gna_ttxr*m_ttxr*m_ttxr*m_ttxr*hsr*nao*F*((exp((v-enattxr)*F/(R*T))-1.0));

/**If**************************************************************************/

y_inf = 1.0/(1.0 + exp((v+vhalf_gh)/16.3));
tau_y_1_2 = 1.5049/(exp(-(v+590.3)*0.01094)+ exp((v-85.1)/17.2));
y_1_2 = y_1_2 + ddt*(y_inf - y_1_2)/tau_y_1_2;
ihk  = 0.6167*gh*y_1_2*(v - ek);
ihna = 0.3833*gh*y_1_2*(v - ena);
ih = (ihk + ihna);

/*Ito*************************************************************************/
q_inf = 1.0/(1.0+exp((v+49.0)/13.0));
tau_q = (6.06 + 39.102/(0.57*exp(-0.08*(v+44.0))+0.065*exp(0.1*(v+45.93))))/0.67; 
q = q + ddt*((q_inf-q)/tau_q);
r_inf = 1.0/(1.0+exp(-(v-19.3)/15.0));
tau_r = (2.75+14.40516/(1.037*exp(0.09*(v+30.61))+0.369*exp(-0.12*(v+23.84))))/0.303;
r = r + ddt*((r_inf-r)/tau_r);
ito = gto*q*r*(v-ek);
 
/*Isus***********************************************************************/
isus = gsus*r*(v-ek);
/*Inak***********************************************************************/
inak = inakmax*(pow(ko,1.2)/(pow(kmkp,1.2)+pow(ko,1.2)))*(pow(nai,1.3)/(pow(kmnap,1.3)+pow(nai,1.3)))/(1.0+exp(-(v-ena+120.0)/30.0));

/****iNaCa*******************************************************************/

di=1+(casub/Kci)*(1+exp(-Qci*v*FRT)+nai/Kcni)+(nai/K1ni)*(1+(nai/K2ni)*(1+nai/K3ni));
doo=1+(cao/Kco)*(1+exp(Qco*v*FRT))+(nao/K1no)*(1+(nao/K2no)*(1+nao/K3no));
k43=nai/(K3ni+nai);
k12=(casub/Kci)*exp(-Qci*v*FRT)/di;
k14=(nai/K1ni)*(nai/K2ni)*(1+nai/K3ni)*exp(Qn*v*FRT/2.0)/di;
k41=exp(-Qn*v*FRT/2.0);
k34=nao/(K3no+nao);
k21=(cao/Kco)*exp(Qco*v*FRT)/doo;
k23=(nao/K1no)*(nao/K2no)*(1+nao/K3no)*exp(-Qn*v*FRT/2.0)/doo;
k32=exp(Qn*v*FRT/2);
x1=k34*k41*(k23+k21)+k21*k32*(k43+k41);
x2=k43*k32*(k14+k12)+k41*k12*(k34+k32);
x3=k43*k14*(k23+k21)+k12*k23*(k43+k41);
x4=k34*k23*(k14+k12)+k21*k14*(k34+k32);
inaca = kNaCa*(k21*x2-k12*x1)/(x1+x2+x3+x4);   
ca_flux = (ical12+ical13+icat-2.0*inaca+ibca)/(2.0*F);
Jcadif = (casub - cai)/tdifca;
kcasr = maxsr - (maxsr - minsr)/(1.0 + pow(eca50sr/carel,hsrr));
kosrca = koca/kcasr;
kisrca = kica*kcasr;
resting = resting + ddt*(kim*resting_inactivated - kisrca*casub*resting - kosrca*casub*casub*resting + kom*open);
open    = open    + ddt*(kosrca*casub*casub*resting - kom*open - kisrca*casub*open + kim*inactivated);
inactivated = inactivated + ddt*(kisrca*casub*open - kim*inactivated - kom*inactivated + kosrca*casub*casub*resting_inactivated);
resting_inactivated = resting_inactivated + ddt*(kom*inactivated - kosrca*casub*casub*resting_inactivated - kim*resting_inactivated + kisrca*casub*resting);
Jrel = ks*open*(carel - casub);
Jup = Pup*(pow(cai/pumpkmf,pumphill) - pow(caup/pumpkmr,pumphill))/(1.0 + pow(cai/pumpkmf,pumphill) + pow(caup/pumpkmr,pumphill));
Jtr  = (caup - carel)/Ttr;
dFtc  = kfTC*cai*(1.0-Ftc)-kbTC*Ftc;
dFtmc = kfTMC*cai*(1.0-Ftmc-Ftmm)-kbTMC*Ftmc;
dFtmm = kfTMM*Mgi*(1.0-Ftmc-Ftmm)-kbTMM*Ftmm;
dFcms = kfCM*casub*(1.0-Fcms)-kbCM*Fcms;
dFcmi = kfCM*cai*(1.0-Fcmi)-kbCM*Fcmi;
dFcq  = kfCQ*carel*(1.0-Fcq)-kbCQ*Fcq;
Ftc = Ftc + ddt*dFtc;        
Ftmc = Ftmc + ddt*dFtmc;     
Ftmm = Ftmm + ddt*dFtmm;     
Fcms = Fcms + ddt*dFcms;     
Fcmi = Fcmi + ddt*dFcmi;     
Fcq = Fcq + ddt*dFcq;        
casub = casub + ddt*((-ca_flux+Jrel*vrel)/vsub-Jcadif-ConcCM*dFcms);
cai = cai + ddt*((Jcadif*vsub-Jup*vup)/vi - (ConcCM*dFcmi + ConcTC*dFtc + ConcTMC*dFtmc)); 
carel = carel + ddt*(Jtr - Jrel - ConcCQ*dFcq);
caup = caup + ddt*(Jup-Jtr*vrel/vup);
dvdtold = dvdt;
total_current = ih+ina_ttxr+ina_ttxs+ical12+ical13+iks+ikr+ik1+ist+ib+icat+inak+isus+inaca+ito;
dvdt = - total_current/capacitance;
vnew = v  + ddt*dvdt;
ena = (R*T/F)*log(nao/nai);
ek  = (R*T/F)*log(ko/ki);
eks = ((R*T)/F)*log((ko+0.12*nao)/(ki+0.12*nai));
eca = (R*T/(2*F))*log(cao/casub);
nai_tot = ihna+ina_ttxr+ina_ttxs+3.0*inak+3.0*inaca+ist+ibna;
ki_tot = ihk+iks+ikr+ik1+ibk-2.0*inak+isus+ito;
nai = nai - ddt*(nai_tot)/(F*vi);
ki = ki - ddt*(ki_tot)/(F*vi);

if(dvdt>=0.0&&dvdtold<0.0){
 min_potential[param_counter] = v;
 nai_min = nai;
 ki_min = ki;
 tmin_potential[param_counter] = sstime;
 start_output = 1;
}
if(dvdt>dvdtmax[param_counter]&&start_output>0){
 dvdtmax[param_counter] = dvdt;
 apd_start[param_counter] = sstime;
 vdvdtmax[param_counter] = v;	
}
if(dvdtold>0.0&&dvdt<=0.0){
 max_potential[param_counter] = v;
 top_slope[param_counter] = (max_potential[param_counter]-min_potential[param_counter])/(sstime - tmin_potential[param_counter]);
}
if((param_counter>0)&&(dvdtold<=top_slope[param_counter-1])&&(dvdt>top_slope[param_counter-1])){
 top[param_counter] = v;
 ddr[param_counter] = (v - min_potential[param_counter])/(sstime - tmin_potential[param_counter]);
}
if(vnew<=0.5*min_potential[param_counter]&&v>0.5*min_potential[param_counter]){
 if(apd_start[param_counter]>0.0)
 apd50[param_counter] = sstime - apd_start[param_counter];
}
if(vnew<=0.9*min_potential[param_counter]&&v>0.9*min_potential[param_counter]){
 if(apd_start[param_counter]>0.0){
  apd_end[param_counter] = sstime;
  apd90[param_counter] = sstime - apd_start[param_counter];
  cycle_length[param_counter] = apd_start[param_counter]-apd_start[param_counter-1];
  printf("%10.5f\t",min_potential[param_counter]);
  printf("%10.5f\t",max_potential[param_counter]);
  printf("%10.5f\t",dvdtmax[param_counter]);
  printf("%10.5f\t",apd50[param_counter]);
  printf("%10.5f\t",apd90[param_counter]);
  printf("%10.5f\t",cycle_length[param_counter]);
  printf("%10.5f\t",ddr[param_counter]);
  printf("%10.5f\t",top[param_counter]);
  printf("\n");
  param_counter++;
 }
}
v = vnew;

if(current_trace_counter%10==0){
fprintf(outputcurrents,"%10.10f\t%10.10f\t",sstime,v);
fprintf(outputcurrents,"\n");
}

current_trace_counter++;

} // end of time loop

fclose(outputcurrents);

return 0;
}

Loading data, please wait...