/* * A Model of Sinoatrial Node Vagal Control * * Model Status * * This model is valid CellML. It has been unit checked and curated * and is known to reproduce the published results in COR and OpenCell. * * Model Structure * * Sinoatrial (SA) node cells have an inherent ability to generate * a depolarising, unstable resting potential leading to automaticity. * The rhythmic, electrical activity of the sinoatrial cells set * the rate at which the entire heart beats, hence the sinoatrial * node myocytes are referred to as the pacemaker cells. This cardiac * pacemaker activity is under vagal control. The ionic mechanisms * underlying the vagal inhibition of the cardiac pacemaker are * the subject of investigation in a mathematical model published * by Socrates Dokos, Branko Celler and Nigel Lovell (1996). * * In this paper, the authors review the existing knowledge surrounding * the vagal control of sinoatrial rhythm. It is known that following * vagal stimulation, acetylcholine (ACh) is released into the * parasympathetic neuroeffector junction, and then binds to muscarinic * receptors on the plasma membrane of the SA node cells. This * receptor-binding triggers membrane hyperpolarisation, and/or * decreases the rate of pacemaker depolarisation, in turn prolonging * the spontaneous cycle duration, and decreasing the rate of autorhythmic * firing. The principal mechanism underlying this ACh-mediated * inhibition of the cardiac pacemaker is an increase in the membrane * permeability to K+. However, it has been suggested that the * influence of ACh on other ion currents may also have a significant * effect on pacemaker activity. * * The focus of the Dokos et al. 1996 model was to gain a better * understanding of the mechanisms underlying vagal control of * the cardiac pacemaker. Their model was based on a wide range * of electrophysiological data, and their aim was to reproduce * these experimental results with their mathematical model. This * model is an extension of their previously published mathematical * model of the SA node. In this new model, the background potassium * current ib,K has been replaced by an ACh-activated potassium * current iK,ACh. In addition, the new model incorporates the * influence of ACh on the other ionic currents, such as the inhibition * of the hyperpolarisation-activated current if, and the inhibition * of the L-type calcium current iCa,L. * * Vagal control of pacemaker activity is modelled using a three * compartment model describing ACh release and uptake in the neuroeffector * junction (see below). Upon vagal stimulus, ACh is released into * the neuroeffector junction, activating iK,ACh and inhibiting * if and iCa,L. * * The complete original paper reference is cited below: * * Vagal Control of Sinoatrial Rhythm: a Mathematical Model, Socrates * Dokos, Branko Celler, and Nigel Lovell, 1996, Journal of Theoretical * Biology , 182, 21-44. PubMed ID: 8917735 * * cell diagram * * [[Image file: dokos_ii_1996.png]] * * A schematic diagram of the Dokos et al. 1996 mathematical model * of vagal control of cardiac pacemaker activity. */ import nsrunit; // Warning: unit conversion turned off due to unit errors in 39 equation(s) unit conversion off; unit flux=1 meter^(-3)*second^(-1)*mole^1; unit first_order_rate_constant=1 second^(-1); unit second_order_rate_constant=1 meter^3*second^(-1)*mole^(-1); // unit millivolt predefined unit per_millivolt=1E3 kilogram^(-1)*meter^(-2)*second^3*ampere^1; unit per_millivolt_second=1E3 kilogram^(-1)*meter^(-2)*second^2*ampere^1; unit nanoS=1E-9 kilogram^(-1)*meter^(-2)*second^3*ampere^2; unit picoF=1E-12 kilogram^(-1)*meter^(-2)*second^4*ampere^2; unit picoA_per_millimolar=1E-12 meter^4.77*ampere^1*mole^(-1.59); unit picoA=1E-12 ampere^1; // unit millimolar predefined unit joule_per_mole_kelvin=1 kilogram^1*meter^2*second^(-2)*kelvin^(-1)*mole^(-1); unit coulomb_per_millimole=.001 second^(-1)*ampere^(-1)*mole^1; // unit picolitre predefined math main { //Warning: the following variables were set 'extern' or given // an initial value of '0' because the model would otherwise be // underdetermined: ACh, w, v, u, ACh_ms, ACh_ex realDomain time second; time.min=0; extern time.max; extern time.delta; real R joule_per_mole_kelvin; R=8.32; real T kelvin; T=310.0; real F coulomb_per_millimole; F=96.49; real E(time) millivolt; when(time=time.min) E=-64.9; real C picoF; C=32.0; real i_tot(time) picoA; real i_CaL(time) picoA; real i_CaT(time) picoA; real i_Na(time) picoA; real i_K(time) picoA; real i_f(time) picoA; real i_p(time) picoA; real i_NaCa(time) picoA; real i_bNa(time) picoA; real i_KACh(time) picoA; real E_Ca(time) millivolt; real E_Na(time) millivolt; real E_K(time) millivolt; real Cai(time) millimolar; when(time=time.min) Cai=0.000034; real Cao(time) millimolar; when(time=time.min) Cao=2.0004; real Nai(time) millimolar; when(time=time.min) Nai=7.4994; real Nao(time) millimolar; when(time=time.min) Nao=139.9929; real Ki(time) millimolar; when(time=time.min) Ki=140.0073; real Ko(time) millimolar; when(time=time.min) Ko=5.4243; real g_CaL nanoS; g_CaL=400.0; real ACh(time) millimolar; //Warning: Assuming zero initial condition; nothing provided in original CellML model. when(time=time.min) ACh=0; real L_type_calcium_current.d(time) dimensionless; when(time=time.min) L_type_calcium_current.d=0.0001; real L_type_calcium_current.f(time) dimensionless; when(time=time.min) L_type_calcium_current.f=0.1505; real f2(time) dimensionless; when(time=time.min) f2=0.2190; real w(time) dimensionless; //Warning: Assuming zero initial condition; nothing provided in original CellML model. when(time=time.min) w=0; real L_type_calcium_current_d_gate.d_infinity(time) dimensionless; real L_type_calcium_current_d_gate.tau_d second; L_type_calcium_current_d_gate.tau_d=0.002; real L_type_calcium_current_f_gate.f_infinity(time) dimensionless; real L_type_calcium_current_f_gate.tau_f(time) second; real alpha_f2 first_order_rate_constant; alpha_f2=3.0; real beta_f2 second_order_rate_constant; beta_f2=40000.0; real alpha_w(time) first_order_rate_constant; real beta_w(time) first_order_rate_constant; real g_CaT nanoS; g_CaT=85.0; real T_type_calcium_current.d(time) dimensionless; when(time=time.min) T_type_calcium_current.d=0.0010; real T_type_calcium_current.f(time) dimensionless; when(time=time.min) T_type_calcium_current.f=0.1328; real T_type_calcium_current_d_gate.d_infinity(time) dimensionless; real T_type_calcium_current_d_gate.tau_d(time) second; real T_type_calcium_current_f_gate.f_infinity(time) dimensionless; real T_type_calcium_current_f_gate.tau_f(time) second; real g_Na nanoS; g_Na=250.0; real m(time) dimensionless; when(time=time.min) m=0.0139; real h(time) dimensionless; when(time=time.min) h=0.0087; real alpha_m(time) first_order_rate_constant; real beta_m(time) first_order_rate_constant; real alpha_h(time) first_order_rate_constant; real beta_h(time) first_order_rate_constant; real i_KK(time) picoA; real i_KNa(time) picoA; real Kk picoA_per_millimolar; Kk=0.26; real P_KNa dimensionless; P_KNa=0.035; real x(time) dimensionless; when(time=time.min) x=0.5682; real x_infinity(time) dimensionless; real tau_x(time) second; real i_fNa(time) picoA; real i_fK(time) picoA; real Kmf millimolar; Kmf=10.3; real g_fNa nanoS; g_fNa=8.1; real g_fK nanoS; g_fK=13.5; real y(time) dimensionless; when(time=time.min) y=0.0287; real alpha_y(time) first_order_rate_constant; real beta_y(time) first_order_rate_constant; real v(time) dimensionless; //Warning: Assuming zero initial condition; nothing provided in original CellML model. when(time=time.min) v=0; real alpha_v(time) first_order_rate_constant; real beta_v(time) first_order_rate_constant; real E_shift_max millivolt; E_shift_max=10.0; real KmNa millimolar; KmNa=40.0; real KmK millimolar; KmK=1.0; real i_pmax picoA; i_pmax=226.0; real kNaCa picoA; kNaCa=4000.0; real x1(time) dimensionless; real x2(time) dimensionless; real x3(time) dimensionless; real x4(time) dimensionless; real k41(time) dimensionless; real k34(time) dimensionless; real k23(time) dimensionless; real k21(time) dimensionless; real k32(time) dimensionless; real k43(time) 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.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.0; real K2no millimolar; K2no=561.4; real Kco millimolar; Kco=3.663; real do(time) dimensionless; real di(time) dimensionless; real g_Nab nanoS; g_Nab=0.24; real KbK picoA_per_millimolar; KbK=0.07; real KKACh picoA_per_millimolar; KKACh=2.0; real u(time) dimensionless; //Warning: Assuming zero initial condition; nothing provided in original CellML model. when(time=time.min) u=0; real alpha_u(time) first_order_rate_constant; real beta_u(time) first_order_rate_constant; real ACh_ms(time) millimolar; //Warning: Assuming zero initial condition; nothing provided in original CellML model. when(time=time.min) ACh_ms=0; real ACh_ex(time) millimolar; //Warning: Assuming zero initial condition; nothing provided in original CellML model. when(time=time.min) ACh_ex=0; real k1 dimensionless; k1=0.01; real k2 dimensionless; k2=0.0005; real F_ms dimensionless; F_ms=1.0; real F_ex dimensionless; F_ex=1.0; real k_D first_order_rate_constant; k_D=0.5; real k_E first_order_rate_constant; k_E=0.5; real k_H first_order_rate_constant; k_H=14.0; real i_up(time) picoA; real i_tr(time) picoA; real i_rel(time) picoA; real V_rel picolitre; V_rel=0.015; real V_up picolitre; V_up=0.035; real i_up_max picoA; i_up_max=21.2; real alpha_rel first_order_rate_constant; real alpha_tr first_order_rate_constant; real KmCaup millimolar; KmCaup=0.0005; real KmCarel millimolar; KmCarel=0.001; real tau_rel second; tau_rel=0.005; real tau_tr second; tau_tr=0.4; real Caup(time) millimolar; when(time=time.min) Caup=0.5832; real Carel(time) millimolar; when(time=time.min) Carel=0.1101; real V_i picolitre; V_i=2.5; real V_e picolitre; V_e=0.5; real tau_b second; tau_b=0.1; real Nab millimolar; Nab=140.0; real Cab millimolar; Cab=2.0; real Kb millimolar; Kb=5.4; // // // E:time=((-1)*i_tot/C); i_tot=(i_CaL+i_CaT+i_Na+i_K+i_f+i_p+i_NaCa+i_bNa+i_KACh); // E_Ca=(R*T/(2*F)*ln(Cao/Cai)); E_Na=(R*T/F*ln(Nao/Nai)); E_K=(R*T/F*ln(Ko/Ki)); // i_CaL=(w*g_CaL*L_type_calcium_current.d*L_type_calcium_current.f*f2*(E+(-1)*E_Ca+(75 millivolt))); // L_type_calcium_current_d_gate.d_infinity=(1/(1+exp((E+(6.6 millivolt))/((-6.6 millivolt))))); L_type_calcium_current.d:time=((L_type_calcium_current_d_gate.d_infinity-L_type_calcium_current.d)/L_type_calcium_current_d_gate.tau_d); // L_type_calcium_current_f_gate.f_infinity=(1/(1+exp((E+(25 millivolt))/(6 millivolt)))); L_type_calcium_current_f_gate.tau_f=((.031 second)+(1 second)/(1+exp((E+(37.6 millivolt))/(8.1 millivolt)))); L_type_calcium_current.f:time=((L_type_calcium_current_f_gate.f_infinity-L_type_calcium_current.f)/L_type_calcium_current_f_gate.tau_f); // f2:time=(alpha_f2*(1-f2)-beta_f2*Cai*f2); // alpha_w=(((.0033 first_order_rate_constant)*ACh+6E-7)/(ACh+3.6E-5)); beta_w=((.0037 first_order_rate_constant)*ACh/(ACh+3.6E-5)); w:time=(alpha_w*(1-w)-beta_w*w); // i_CaT=(g_CaT*T_type_calcium_current.d*T_type_calcium_current.f*(E+(-1)*E_Ca+(75 millivolt))); // T_type_calcium_current_d_gate.d_infinity=(1/(1+exp((E+(23 millivolt))/((-6.1 millivolt))))); T_type_calcium_current_d_gate.tau_d=((6E-4 second)+(.0054 second)/(1+exp(.03*(E+(100 millivolt))))); T_type_calcium_current.d:time=((T_type_calcium_current_d_gate.d_infinity-T_type_calcium_current.d)/T_type_calcium_current_d_gate.tau_d); // T_type_calcium_current_f_gate.f_infinity=(1/(1+exp((E+(75 millivolt))/(6.6 millivolt)))); T_type_calcium_current_f_gate.tau_f=((.001 second)+(.04 second)/(1+exp(.08*(E+(65 millivolt))))); T_type_calcium_current.f:time=((T_type_calcium_current_f_gate.f_infinity-T_type_calcium_current.f)/T_type_calcium_current_f_gate.tau_f); // i_Na=(g_Na*m^3*h*(E-E_Na)); // alpha_m=((200 per_millivolt_second)*(E+(34.3 millivolt))/(1-exp((-0.09)*(E+(34.3 millivolt))))); beta_m=((8E3 first_order_rate_constant)*exp((-0.15)*(E+(56.2 millivolt)))); m:time=(alpha_m*(1-m)-beta_m*m); // alpha_h=((32.4 first_order_rate_constant)*exp((-0.14)*(E+(93.4 millivolt)))); beta_h=((709 first_order_rate_constant)/(1+4.2*exp((-0.06)*(E+(45.4 millivolt))))); h:time=(alpha_h*(1-h)-beta_h*h); // i_K=(i_KK+i_KNa); i_KK=(x*Kk*Ko^.59*(Ki-Ko*exp((-1)*E*F/(R*T)))); i_KNa=(x*Kk*P_KNa*Ko^.59*(Nai-Nao*exp((-1)*E*F/(R*T)))); // x_infinity=(1/(1+exp((E+(25.1 millivolt))/((-7.4 millivolt))))); tau_x=((1 second)/(17*exp(.0398*E)+.211*exp((-0.051)*E))); x:time=((x_infinity-x)/tau_x); // i_f=(i_fK+i_fNa); i_fNa=(y*(Ko^1.83/(Ko^1.83+Kmf^1.83))*(g_fNa*(Ko*(E-E_Na)))); i_fK=(y*(Ko^1.83/(Ko^1.83+Kmf^1.83))*(g_fK*(Ko*(E-E_K)))); // alpha_v=((4 first_order_rate_constant)*ACh/(ACh+(8E-5 millimolar))); beta_v=((1E-5 flux)/(ACh+(8E-5 millimolar))); v:time=(alpha_v*(1-v)-beta_v*v); alpha_y=((.36 per_millivolt_second)*(E+(137.8 millivolt)+v^4*E_shift_max)/(exp(.0666*(E+(137.8 millivolt)+v^4*E_shift_max))-1)); beta_y=((.1 per_millivolt_second)*(E+(76.3 millivolt)+v^4*E_shift_max)/(1-exp((-0.21)*(E+(76.3 millivolt)+v^4*E_shift_max)))); y:time=(alpha_y*(1-y)-beta_y*y); // i_p=(i_pmax*(Nai/(Nai+KmNa))*(Ko/(Ko+KmK))*(1-((E-(40 millivolt))/(211 millivolt))^2)); // i_NaCa=(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)); k12=(Cai/Kci*exp((-1)*Qci*E*F/(R*T))/di); k14=((Nai^2/(K1ni*K2ni)+Nai^3/(K1ni*K2ni*K3ni))*exp(Qn*E*F/(R*T))/di); k41=exp((-1)*(Qn*E*F)/(2*R*T)); di=(1+Cai/Kci+Cai/Kci*exp((-1)*(Qci*E*F)/(R*T))+Cai*Nai/(Kci*Kcni)+Nai/K1ni+Nai^2/(K1ni*K2ni)+Nai^3/(K1ni*K2ni*K3ni)); k34=(Nao/(K3no+Nao)); k21=(Cao/Kco*exp(Qco*E*F/(R*T))/do); k23=(Nao^2/(K1no*K2no)*(Nao^3/(K1no*K2no*K3no))*exp((-1)*(Qn*E*F)/(2*R*T))/do); k32=exp(Qn*E*F/(2*R*T)); do=(1+Cao/Kco+Cao/Kco*exp((-1)*(Qco*E*F)/(R*T))+Nao/K1no+Nao^2/(K1no*K2no)+Nao^3/(K1no*K2no*K3no)); // i_bNa=(g_Nab*(E-E_Na)); // i_KACh=((KbK+u^4*KKACh)*Ko^.41*(Ki-Ko*exp((-1)*E*F/(R*T)))); // alpha_u=((4.7 first_order_rate_constant)*ACh/(ACh+(1.2E-4 millimolar))); beta_u=((2.7E-4 flux)/(ACh+(1.2E-4 millimolar))); u:time=(alpha_u*(1-u)-beta_u*u); // ACh_ms:time=(k_H/F_ms*ACh+k_E/F_ms*ACh_ex); ACh:time=((-1)*k_H*ACh-k_D*(ACh-ACh_ex)); ACh_ex:time=((-1)*(k_E/F_ex)*ACh-k_D/F_ms*(ACh_ex-ACh)); // i_up=(i_up_max*(Cai^2/(Cai^2+KmCaup^2))); i_tr=(alpha_tr*Caup); i_rel=(alpha_rel*Carel*(Cai^2/(Cai^2+KmCarel^2))); alpha_rel=(2*V_rel*F/tau_rel); alpha_tr=(2*V_rel*F/tau_rel); // Nai:time=((-1)*(i_bNa+i_fNa+i_Na+3*i_p+3*i_NaCa+i_KNa)/(F*V_i)); Nao:time=((i_bNa+i_fNa+i_Na+3*i_p+3*i_NaCa+i_KNa)/(F*V_e)+(Nab-Nao)/tau_b); Ki:time=((-1)*(i_KK+i_fK+(-2)*i_p+i_KACh)/(F*V_i)); Ko:time=((i_KK+i_fK+(-2)*i_p+i_KACh)/(F*V_e)+(Kb-Ko)/tau_b); Cai:time=((-1)*(i_CaL+i_CaT+(-2)*i_NaCa+i_up+(-1)*i_rel)/(2*F*V_i)); Cao:time=((i_CaL+i_CaT+(-2)*i_NaCa)/(2*F*V_e)+(Cab-Cao)/tau_b); Caup:time=((i_up-i_tr)/(2*V_up*F)); Carel:time=((i_tr-i_rel)/(2*V_rel*F)); }