/* * Synergism of coupled subsarcolemmal Ca2+ clocks and sarcolemmal * voltage clocks confers robust and flexible pacemaker function * in a novel pacemaker cell model * * Model Status * * This CellML model has been tested in OpenCell and COR and the * units are balanced. The model runs to recreate the published * results. This particular variant of the model uses the initial * conditions from the author's original code to create a steady * state model. * * Model Structure * * ABSTRACT: Recent experimental studies have demonstrated that * sinoatrial node cells (SANC) generate spontaneous, rhythmic, * local subsarcolemmal Ca(2+) releases (Ca(2+) clock), which occur * during late diastolic depolarization (DD) and interact with * the classic sarcolemmal voltage oscillator (membrane clock) * by activating Na(+)-Ca(2+) exchanger current (I(NCX)). This * and other interactions between clocks, however, are not captured * by existing essentially membrane-delimited cardiac pacemaker * cell numerical models. Using wide-scale parametric analysis * of classic formulations of membrane clock and Ca(2+) cycling, * we have constructed and initially explored a prototype rabbit * SANC model featuring both clocks. Our coupled oscillator system * exhibits greater robustness and flexibility than membrane clock * operating alone. Rhythmic spontaneous Ca(2+) releases of sarcoplasmic * reticulum (SR)-based Ca(2+) clock ignite rhythmic action potentials * via late DD I(NCX) over much broader ranges of membrane clock * parameters [e.g., L-type Ca(2+) current (I(CaL)) and/or hyperpolarization-activated * ("funny") current (I(f)) conductances]. The system Ca(2+) clock * includes SR and sarcolemmal Ca(2+) fluxes, which optimize cell * Ca(2+) balance to increase amplitudes of both SR Ca(2+) release * and late DD I(NCX) as SR Ca(2+) pumping rate increases, resulting * in a broad pacemaker rate modulation (1.8-4.6 Hz). In contrast, * the rate modulation range via membrane clock parameters is substantially * smaller when Ca(2+) clock is unchanged or lacking. When Ca(2+) * clock is disabled, the system parametric space for fail-safe * SANC operation considerably shrinks: without rhythmic late DD * I(NCX) ignition signals membrane clock substantially slows, * becomes dysrhythmic, or halts. In conclusion, the Ca(2+) clock * is a new critical dimension in SANC function. A synergism of * the coupled function of Ca(2+) and membrane clocks confers fail-safe * SANC operation at greatly varying rates. * * The original paper reference is cited below: * * Synergism of coupled subsarcolemmal Ca2+ clocks and sarcolemmal * voltage clocks confers robust and flexible pacemaker function * in a novel pacemaker cell model, Victor A. Maltsev and Edward * G. Lakatta, 2009, American Journal of Physiology, 286, H594-H615. * PubMed ID: 19136600 * * reaction diagram * * [[Image file: maltsev_2009.png]] * * Schematic diagram of the interacting Ca2+ clock and membrane * clock in a model of rabbit sinoatrial node cells. */ import nsrunit; unit conversion on; unit per_millisecond=1E3 second^(-1); unit millimolar_per_millisecond=1E3 meter^(-3)*second^(-1)*mole^1; unit per_millimolar_millisecond=1E3 meter^3*second^(-1)*mole^(-1); unit per_millimolar2_millisecond=1E3 meter^6*second^(-1)*mole^(-2); // unit millisecond predefined // unit millivolt predefined unit per_millivolt=1E3 kilogram^(-1)*meter^(-2)*second^3*ampere^1; unit per_millivolt_millisecond=1E6 kilogram^(-1)*meter^(-2)*second^2*ampere^1; unit nanoS_per_picoF=1E3 second^(-1); unit picoF=1E-12 kilogram^(-1)*meter^(-2)*second^4*ampere^2; unit per_picoF=1E12 kilogram^1*meter^2*second^(-4)*ampere^(-2); unit picoA_per_picoF=1 kilogram^1*meter^2*second^(-4)*ampere^(-1); unit picoA_per_millimolar_picoF=1 kilogram^1*meter^5*second^(-4)*ampere^(-1)*mole^(-1); unit picoA=1E-12 ampere^1; // unit picolitre predefined unit picolitre_per_femtolitre=1E3 dimensionless; // unit micrometre predefined // unit millimolar predefined unit joule_per_kilomole_kelvin=.001 kilogram^1*meter^2*second^(-2)*kelvin^(-1)*mole^(-1); unit coulomb_per_mole=1 second^1*ampere^1*mole^(-1); math main { realDomain time millisecond; time.min=0; extern time.max; extern time.delta; real Vm(time) millivolt; when(time=time.min) Vm=-57.9639346865; real Cm picoF; Cm=32; real i_CaT(time) picoA; real i_CaL(time) picoA; real i_f(time) picoA; real i_st(time) picoA; real i_Kr(time) picoA; real i_Ks(time) picoA; real i_to(time) picoA; real i_sus(time) picoA; real i_NaK(time) picoA; real i_NaCa(time) picoA; real i_b_Ca(time) picoA; real i_b_Na(time) picoA; real E_Na millivolt; real E_K millivolt; real E_Ks millivolt; real electric_potentials.R joule_per_kilomole_kelvin; electric_potentials.R=8314.4; real T kelvin; T=310.15; real F coulomb_per_mole; F=96485; real Nai millimolar; Nai=10; real Nao millimolar; Nao=140; real Ki millimolar; Ki=140; real Ko millimolar; Ko=5.4; real E_CaL millivolt; E_CaL=45; real g_CaL nanoS_per_picoF; g_CaL=0.464; real Ca_sub(time) millimolar; when(time=time.min) Ca_sub=0.000138112560112; real dL(time) dimensionless; when(time=time.min) dL=0.000584545564405; real fL(time) dimensionless; when(time=time.min) fL=0.862381249774; real fCa(time) dimensionless; when(time=time.min) fCa=0.711395919653; real dL_infinity(time) dimensionless; real tau_dL(time) millisecond; real alpha_dL(time) per_millisecond; real beta_dL(time) per_millisecond; real adVm(time) millivolt; real bdVm(time) millivolt; real fL_infinity(time) dimensionless; real tau_fL(time) millisecond; real alpha_fCa per_millisecond; alpha_fCa=0.021; real fCa_infinity(time) dimensionless; real tau_fCa(time) millisecond; real Km_fCa millimolar; Km_fCa=0.00035; real g_CaT nanoS_per_picoF; g_CaT=0.1832; real E_CaT millivolt; E_CaT=45; real dT(time) dimensionless; when(time=time.min) dT=0.00504393374639; real fT(time) dimensionless; when(time=time.min) fT=0.420757825415; real dT_infinity(time) dimensionless; real tau_dT(time) millisecond; real fT_infinity(time) dimensionless; real tau_fT(time) millisecond; real g_Kr nanoS_per_picoF; g_Kr=0.08113973; real paS(time) dimensionless; when(time=time.min) paS=0.453100576739; real paF(time) dimensionless; when(time=time.min) paF=0.144755091176; real pi_(time) dimensionless; when(time=time.min) pi_=0.849409822329; real pa_infinity(time) dimensionless; real tau_paS(time) millisecond; real tau_paF(time) millisecond; real pi_infinity(time) dimensionless; real tau_pi(time) millisecond; real g_Ks nanoS_per_picoF; g_Ks=0.0259; real n(time) dimensionless; when(time=time.min) n=0.0264600410928; real n_infinity(time) dimensionless; real tau_n(time) millisecond; real alpha_n(time) per_millisecond; real beta_n(time) per_millisecond; real g_to nanoS_per_picoF; g_to=0.252; real g_sus nanoS_per_picoF; g_sus=0.02; real q(time) dimensionless; when(time=time.min) q=0.694241313965; real r(time) dimensionless; when(time=time.min) r=0.00558131733359; real q_infinity(time) dimensionless; real tau_q(time) millisecond; real r_infinity(time) dimensionless; real tau_r(time) millisecond; real i_f_Na(time) picoA; real i_f_K(time) picoA; real g_if nanoS_per_picoF; g_if=0.15; real y(time) dimensionless; when(time=time.min) y=0.113643187247; real y_infinity(time) dimensionless; real tau_y(time) millisecond; real VIf_half millivolt; VIf_half=-64; real g_st nanoS_per_picoF; g_st=0.003; real E_st millivolt; E_st=37.4; real qa(time) dimensionless; when(time=time.min) qa=0.42380243163; real qi(time) dimensionless; when(time=time.min) qi=0.447294008304; real qa_infinity(time) dimensionless; real tau_qa(time) millisecond; real alpha_qa(time) per_millisecond; real beta_qa(time) per_millisecond; real qi_infinity(time) dimensionless; real tau_qi(time) millisecond; real alpha_qi(time) per_millisecond; real beta_qi(time) per_millisecond; real g_b_Na nanoS_per_picoF; g_b_Na=0.00486; real Km_Kp millimolar; Km_Kp=1.4; real Km_Nap millimolar; Km_Nap=14; real i_NaK_max picoA_per_picoF; i_NaK_max=2.88; real g_b_Ca nanoS_per_picoF; g_b_Ca=0.0006; real kNaCa picoA_per_picoF; kNaCa=187.5; real x1(time) dimensionless; real x2(time) dimensionless; real x3(time) dimensionless; real x4(time) dimensionless; real k41(time) dimensionless; real k34 dimensionless; real k23(time) dimensionless; real k21(time) dimensionless; real k32(time) dimensionless; real k43 dimensionless; real k12(time) dimensionless; real k14(time) dimensionless; real Qci dimensionless; Qci=0.1369; real Qn dimensionless; Qn=0.4315; real Qco dimensionless; Qco=0; real K3ni millimolar; K3ni=26.44; real Kci millimolar; Kci=0.0207; real K1ni millimolar; K1ni=395.3; real K2ni millimolar; K2ni=2.289; real Kcni millimolar; Kcni=26.44; real K3no millimolar; K3no=4.663; real K1no millimolar; K1no=1628; real K2no millimolar; K2no=561.4; real Kco millimolar; Kco=3.663; real RTOnF millivolt; real do(time) dimensionless; real di(time) dimensionless; real Cao millimolar; Cao=2; real j_SRCarel(time) millimolar_per_millisecond; real j_SRCarel.R(time) dimensionless; when(time=time.min) j_SRCarel.R=0.688047760973; real O(time) dimensionless; when(time=time.min) O=1.7340201253e-7; real I(time) dimensionless; when(time=time.min) I=7.86181717518e-8; real RI(time) dimensionless; when(time=time.min) RI=0.311951987007; real ks per_millisecond; ks=250000; real MaxSR dimensionless; MaxSR=15; real MinSR dimensionless; MinSR=1; real EC50_SR millimolar; EC50_SR=0.45; real HSR dimensionless; HSR=2.5; real koSRCa(time) per_millimolar2_millisecond; real kiSRCa(time) per_millimolar_millisecond; real koCa per_millimolar2_millisecond; koCa=10; real kiCa per_millimolar_millisecond; kiCa=0.5; real kCaSR(time) dimensionless; real kim per_millisecond; kim=0.005; real kom per_millisecond; kom=0.06; real Ca_jsr(time) millimolar; when(time=time.min) Ca_jsr=0.316762674605; real j_Ca_dif(time) millimolar_per_millisecond; real j_up(time) millimolar_per_millisecond; real j_tr(time) millimolar_per_millisecond; real tau_dif_Ca millisecond; tau_dif_Ca=0.04; real tau_tr millisecond; tau_tr=40; real P_up millimolar_per_millisecond; P_up=0.012; real K_up millimolar; K_up=0.0006; real Ca_nsr(time) millimolar; when(time=time.min) Ca_nsr=1.49348117734; real Cai(time) millimolar; when(time=time.min) Cai=0.000150018670943; real TC_tot millimolar; TC_tot=0.031; real TMC_tot millimolar; TMC_tot=0.062; real CM_tot millimolar; CM_tot=0.045; real CQ_tot millimolar; CQ_tot=10; real delta_fTC(time) per_millisecond; real delta_fTMC(time) per_millisecond; real delta_fCMs(time) per_millisecond; real delta_fCMi(time) per_millisecond; real delta_fCQ(time) per_millisecond; real delta_fTMM(time) per_millisecond; real fTMM(time) dimensionless; when(time=time.min) fTMM=0.501049376634; real fCMi(time) dimensionless; when(time=time.min) fCMi=0.0594880901438; real fCMs(time) dimensionless; when(time=time.min) fCMs=0.054381370046; real fTC(time) dimensionless; when(time=time.min) fTC=0.0291316176172; real fTMC(time) dimensionless; when(time=time.min) fTMC=0.432694959597; real fCQ(time) dimensionless; when(time=time.min) fCQ=0.273207128393; real kf_TC per_millimolar_millisecond; kf_TC=88.8; real kf_TMM per_millimolar_millisecond; kf_TMM=2.277; real kf_TMC per_millimolar_millisecond; kf_TMC=227.7; real kf_CM per_millimolar_millisecond; kf_CM=227.7; real kf_CQ per_millimolar_millisecond; kf_CQ=0.534; real kb_TC per_millisecond; kb_TC=0.446; real kb_TMC per_millisecond; kb_TMC=0.00751; real kb_TMM per_millisecond; kb_TMM=0.751; real kb_CM per_millisecond; kb_CM=0.542; real kb_CQ per_millisecond; kb_CQ=0.445; real Mgi millimolar; Mgi=2.5; real V_i picolitre; real V_jsr picolitre; real V_nsr picolitre; real V_sub picolitre; real V_cell picolitre; real V_jsr_part dimensionless; V_jsr_part=0.0012; real V_i_part dimensionless; V_i_part=0.46; real V_nsr_part dimensionless; V_nsr_part=0.0116; real R_cell micrometre; R_cell=4; real L_cell micrometre; L_cell=70; real L_sub micrometre; L_sub=0.02; // // Vm:time=((-1)*(i_CaL+i_CaT+i_f+i_st+i_Kr+i_Ks+i_to+i_sus+i_NaK+i_NaCa+i_b_Ca+i_b_Na)/Cm); // E_Na=(electric_potentials.R*T/F*ln(Nao/Nai)); E_K=(electric_potentials.R*T/F*ln(Ko/Ki)); E_Ks=(electric_potentials.R*T/F*ln((Ko+.12*Nao)/(Ki+.12*Nai))); // i_CaL=(Cm*g_CaL*(Vm-E_CaL)*dL*fL*fCa); // dL_infinity=(1/(1+exp((-1)*(Vm+(13.5 millivolt))/(6 millivolt)))); tau_dL=(1/(alpha_dL+beta_dL)); alpha_dL=((-1)*(.02839 per_millivolt_millisecond)*(adVm+(35 millivolt))/(exp((-1)*(adVm+(35 millivolt))/(2.5 millivolt))-1)-(.0849 per_millivolt_millisecond)*adVm/(exp((-1)*adVm/(4.8 millivolt))-1)); adVm=(if (Vm=((-1)*(35 millivolt))) (-1)*(35.00001 millivolt) else if (Vm=(0 millivolt)) (1E-5 millivolt) else Vm); beta_dL=((.01143 per_millivolt_millisecond)*(bdVm-(5 millivolt))/(exp((bdVm-(5 millivolt))/(2.5 millivolt))-1)); bdVm=(if (Vm=(5 millivolt)) (5.00001 millivolt) else Vm); dL:time=((dL_infinity-dL)/tau_dL); // fL_infinity=(1/(1+exp((Vm+(35 millivolt))/(7.3 millivolt)))); tau_fL=((44.3 millisecond)+(257.1 millisecond)*exp((-1)*((Vm+(32.5 millivolt))/(13.9 millivolt))^2)); fL:time=((fL_infinity-fL)/tau_fL); // fCa_infinity=(Km_fCa/(Km_fCa+Ca_sub)); tau_fCa=(fCa_infinity/alpha_fCa); fCa:time=((fCa_infinity-fCa)/tau_fCa); // i_CaT=(Cm*g_CaT*(Vm-E_CaT)*dT*fT); // dT_infinity=(1/(1+exp((-1)*(Vm+(26.3 millivolt))/(6 millivolt)))); tau_dT=((1 millisecond)/(1.068*exp((Vm+(26.3 millivolt))/(30 millivolt))+1.068*exp((-1)*(Vm+(26.3 millivolt))/(30 millivolt)))); dT:time=((dT_infinity-dT)/tau_dT); // fT_infinity=(1/(1+exp((Vm+(61.7 millivolt))/(5.6 millivolt)))); tau_fT=((1 millisecond)/(.0153*exp((-1)*(Vm+(61.7 millivolt))/(83.3 millivolt))+.015*exp((Vm+(61.7 millivolt))/(15.38 millivolt)))); fT:time=((fT_infinity-fT)/tau_fT); // i_Kr=(Cm*g_Kr*(Vm-E_K)*(.6*paF+.4*paS)*pi_); // pa_infinity=(1/(1+exp((-1)*(Vm+(23.2 millivolt))/(10.6 millivolt)))); tau_paS=((.84655354 millisecond)/(.0042*exp(Vm/(17 millivolt))+1.5E-4*exp((-1)*Vm/(21.6 millivolt)))); tau_paF=((.84655354 millisecond)/(.0372*exp(Vm/(15.9 millivolt))+9.6E-4*exp((-1)*Vm/(22.5 millivolt)))); paS:time=((pa_infinity-paS)/tau_paS); paF:time=((pa_infinity-paF)/tau_paF); // pi_infinity=(1/(1+exp((Vm+(28.6 millivolt))/(17.1 millivolt)))); tau_pi=((1 millisecond)/(.1*exp((-1)*Vm/(54.645 millivolt))+.656*exp(Vm/(106.157 millivolt)))); pi_:time=((pi_infinity-pi_)/tau_pi); // i_Ks=(Cm*g_Ks*(Vm-E_Ks)*n^2); // n_infinity=(alpha_n/(alpha_n+beta_n)); tau_n=(1/(alpha_n+beta_n)); alpha_n=((.014 per_millisecond)/(1+exp((-1)*(Vm-(40 millivolt))/(9 millivolt)))); beta_n=((.001 per_millisecond)*exp((-1)*Vm/(45 millivolt))); n:time=((n_infinity-n)/tau_n); // i_to=(Cm*g_to*(Vm-E_K)*q*r); i_sus=(Cm*g_sus*(Vm-E_K)*r); // q_infinity=(1/(1+exp((Vm+(49 millivolt))/(13 millivolt)))); tau_q=((6.06 millisecond)+(39.102 millisecond)/(.57*exp((-1)*(.08 per_millivolt)*(Vm+(44 millivolt)))+.065*exp((.1 per_millivolt)*(Vm+(45.93 millivolt))))); q:time=((q_infinity-q)/tau_q); // r_infinity=(1/(1+exp((-1)*(Vm-(19.3 millivolt))/(15 millivolt)))); tau_r=((2.75352 millisecond)+(14.40516 millisecond)/(1.037*exp((.09 per_millivolt)*(Vm+(30.61 millivolt)))+.369*exp((-1)*(.12 per_millivolt)*(Vm+(23.84 millivolt))))); r:time=((r_infinity-r)/tau_r); // i_f_Na=(Cm*.3833*g_if*(Vm-E_Na)*y^2); i_f_K=(Cm*.6167*g_if*(Vm-E_K)*y^2); i_f=(i_f_Na+i_f_K); // y_infinity=(1/(1+exp((Vm-VIf_half)/(13.5 millivolt)))); tau_y=((.7166529 millisecond)/(exp((-1)*(Vm+(386.9 millivolt))/(45.302 millivolt))+exp((Vm-(73.08 millivolt))/(19.231 millivolt)))); y:time=((y_infinity-y)/tau_y); // i_st=(Cm*g_st*(Vm-E_st)*qa*qi); // qa_infinity=(1/(1+exp((-1)*(Vm+(57 millivolt))/(5 millivolt)))); tau_qa=(1/(alpha_qa+beta_qa)); alpha_qa=(1/((.15 millisecond)*exp((-1)*Vm/(11 millivolt))+(.2 millisecond)*exp((-1)*Vm/(700 millivolt)))); beta_qa=(1/((16 millisecond)*exp(Vm/(8 millivolt))+(15 millisecond)*exp(Vm/(50 millivolt)))); qa:time=((qa_infinity-qa)/tau_qa); // qi_infinity=(alpha_qi/(alpha_qi+beta_qi)); tau_qi=(6.65/(alpha_qi+beta_qi)); alpha_qi=(1/((3100 millisecond)*exp(Vm/(13 millivolt))+(700 millisecond)*exp(Vm/(70 millivolt)))); beta_qi=(1/((95 millisecond)*exp((-1)*Vm/(10 millivolt))+(50 millisecond)*exp((-1)*Vm/(700 millivolt)))+(2.29E-4 per_millisecond)/(1+exp((-1)*Vm/(5 millivolt)))); qi:time=((qi_infinity-qi)/tau_qi); // i_b_Na=(Cm*g_b_Na*(Vm-E_Na)); // i_NaK=(Cm*i_NaK_max/((1+(Km_Kp/Ko)^1.2)*(1+(Km_Nap/Nai)^1.3)*(1+exp((-1)*(Vm-E_Na+(120 millivolt))/(30 millivolt))))); // i_b_Ca=(Cm*g_b_Ca*(Vm-E_CaL)); // i_NaCa=(Cm*kNaCa*(x2*k21-x1*k12)/(x1+x2+x3+x4)); x1=(k41*k34*(k23+k21)+k21*k32*(k43+k41)); x2=(k32*k43*(k14+k12)+k41*k12*(k34+k32)); x3=(k14*k43*(k23+k21)+k12*k23*(k43+k41)); x4=(k23*k34*(k14+k12)+k14*k21*(k34+k32)); k43=(Nai/(K3ni+Nai)); RTOnF=(electric_potentials.R*T/F); k12=(Ca_sub/Kci*exp((-1)*Qci*Vm/RTOnF)/di); k14=(Nai/K1ni*Nai/K2ni*(1+Nai/K3ni)*exp(Qn*Vm/(2*RTOnF))/di); k41=exp((-1)*Qn*Vm/(2*RTOnF)); di=(1+Ca_sub/Kci*(1+exp((-1)*Qci*Vm/RTOnF)+Nai/Kcni)+Nai/K1ni*(1+Nai/K2ni*(1+Nai/K3ni))); k34=(Nao/(K3no+Nao)); k21=(Cao/Kco*exp(Qco*Vm/RTOnF)/do); k23=(Nao/K1no*Nao/K2no*(1+Nao/K3no)*exp((-1)*Qn*Vm/(2*RTOnF))/do); k32=exp(Qn*Vm/(2*RTOnF)); do=(1+Cao/Kco*(1+exp(Qco*Vm/RTOnF))+Nao/K1no*(1+Nao/K2no*(1+Nao/K3no))); // j_SRCarel=(ks*O*(Ca_jsr-Ca_sub)); kCaSR=(MaxSR-(MaxSR-MinSR)/(1+(EC50_SR/Ca_jsr)^HSR)); koSRCa=(koCa/kCaSR); kiSRCa=(kiCa*kCaSR); j_SRCarel.R:time=(kim*RI-kiSRCa*Ca_sub*j_SRCarel.R-(koSRCa*Ca_sub^2*j_SRCarel.R-kom*O)); O:time=(koSRCa*Ca_sub^2*j_SRCarel.R-kom*O-(kiSRCa*Ca_sub*O-kim*I)); I:time=(kiSRCa*Ca_sub*O-kim*I-(kom*I-koSRCa*Ca_sub^2*RI)); RI:time=(kom*I-koSRCa*Ca_sub^2*RI-(kim*RI-kiSRCa*Ca_sub*j_SRCarel.R)); // j_Ca_dif=((Ca_sub-Cai)/tau_dif_Ca); j_up=(P_up/(1+K_up/Cai)); j_tr=((Ca_nsr-Ca_jsr)/tau_tr); // fTC:time=delta_fTC; delta_fTC=(kf_TC*Cai*(1-fTC)-kb_TC*fTC); fTMC:time=delta_fTMC; delta_fTMC=(kf_TMC*Cai*(1-(fTMC+fTMM))-kb_TMC*fTMC); fTMM:time=delta_fTMM; delta_fTMM=(kf_TMM*Mgi*(1-(fTMC+fTMM))-kb_TMM*fTMM); fCMi:time=delta_fCMi; delta_fCMi=(kf_CM*Cai*(1-fCMi)-kb_CM*fCMi); fCMs:time=delta_fCMs; delta_fCMs=(kf_CM*Ca_sub*(1-fCMs)-kb_CM*fCMs); fCQ:time=delta_fCQ; delta_fCQ=(kf_CQ*Ca_jsr*(1-fCQ)-kb_CQ*fCQ); // Cai:time=((j_Ca_dif*V_sub-j_up*V_nsr)/V_i-(CM_tot*delta_fCMi+TC_tot*delta_fTC+TMC_tot*delta_fTMC)); Ca_sub:time=(j_SRCarel*V_jsr/V_sub-((i_CaL+i_CaT+i_b_Ca-2*i_NaCa)/(2*F*V_sub)+j_Ca_dif+CM_tot*delta_fCMs)); Ca_nsr:time=(j_up-j_tr*V_jsr/V_nsr); Ca_jsr:time=(j_tr-(j_SRCarel+CQ_tot*delta_fCQ)); // V_cell=((.001 picolitre_per_femtolitre)*3.141592653589793*R_cell^2*L_cell); V_sub=((.001 picolitre_per_femtolitre)*2*3.141592653589793*L_sub*(R_cell-L_sub/2)*L_cell); V_jsr=(V_jsr_part*V_cell); V_i=(V_i_part*V_cell-V_sub); V_nsr=(V_nsr_part*V_cell); }