/* * A Simplified Ventricular Myocyte Model * * Model Status * * This CellML model runs in OpenCell and COR to reproduce the * published results. This particular version of the CellML model * has had a stimulus protocol added to allow it to simulate trains * of action potentials. The parameter values used in this variant * (GP) of the Fenton-Karma model are the GP values - the experimental * steady-state curves extracted by Girouard et al. (Girouard et * al., Circulation 93, pp. 603 to 613, 1996) from optical recordings * of membrane potentials on the epicardial surface of the left * ventricle of a guinea pig during plane wave pacing at fixed * cycle lengths - see Table 1 of the 1998 model errata. Simulations * of this CellML model can be run using CMISS. * * Model Structure * * ABSTRACT: Wave propagation in ventricular muscle is rendered * highly anisotropic by the intramural rotation of the fiber. * This rotational anisotropy is especially important because it * can produce a twist of electrical vortices, which measures the * rate of rotation (in degree/mm) of activation wavefronts in * successive planes perpendicular to a line of phase singularity, * or filament. This twist can then significantly alter the dynamics * of the filament. This paper explores this dynamics via numerical * simulation. After a review of the literature, we present modeling * tools that include: (i) a simplified ionic model with three * membrane currents that approximates well the restitution properties * and spiral wave behavior of more complex ionic models of cardiac * action potential (Beeler-Reuter and others), and (ii) a semi-implicit * algorithm for the fast solution of monodomain cable equations * with rotational anisotropy. We then discuss selected results * of a simulation study of vortex dynamics in a parallelepipedal * slab of ventricular muscle of varying wall thickness (S) and * fiber rotation rate (theta(z)). The main finding is that rotational * anisotropy generates a sufficiently large twist to destabilize * a single transmural filament and cause a transition to a wave * turbulent state characterized by a high density of chaotically * moving filaments. This instability is manifested by the propagation * of localized disturbances along the filament and has no previously * known analog in isotropic excitable media. These disturbances * correspond to highly twisted and distorted regions of filament, * or "twistons," that create vortex rings when colliding with * the natural boundaries of the ventricle. Moreover, when sufficiently * twisted, these rings expand and create additional filaments * by further colliding with boundaries. This instability mechanism * is distinct from the commonly invoked patchy failure or wave * breakup that is not observed here during the initial instability. * For modified Beeler-Reuter-like kinetics with stable reentry * in two dimensions, decay into turbulence occurs in the left * ventricle in about one second above a critical wall thickness * in the range of 4-6 mm that matches experiment. However this * decay is suppressed by uniformly decreasing excitability. Specific * experiments to test these results, and a method to characterize * the filament density during fibrillation are discussed. Results * are contrasted with other mechanisms of fibrillation and future * prospects are summarized. * * The original paper reference is cited below: * * Vortex dynamics in three-dimensional continuous myocardium with * fiber rotation: Filament instability and fibrillation, Flavio * Fenton and Alain Karma, 1998, Chaos, 8, 20-47. PubMed ID: 12779708 * * cell diagram * * [[Image file: fenton_1998a.png]] * * A schematic diagram of the three ionic currents described by * the Fenton-Karma model of a ventricular myocyte. */ import nsrunit; unit conversion on; unit ms=.001 second^1; unit per_ms=1E3 second^(-1); unit mV=.001 kilogram^1*meter^2*second^(-3)*ampere^(-1); unit per_mV=1E3 kilogram^(-1)*meter^(-2)*second^3*ampere^1; unit per_mV_ms=1E6 kilogram^(-1)*meter^(-2)*second^2*ampere^1; unit mS_per_cm2=10 kilogram^(-1)*meter^(-4)*second^3*ampere^2; unit uF_per_cm2=.01 kilogram^(-1)*meter^(-4)*second^4*ampere^2; unit uA_per_cm2=.01 meter^(-2)*ampere^1; unit concentration_units=1 meter^(-3)*mole^1; math main { realDomain time ms; time.min=0; extern time.max; extern time.delta; real u(time) dimensionless; when(time=time.min) u=0; real Cm uF_per_cm2; Cm=1; real Vm(time) mV; real V_0 mV; V_0=-85; real V_fi mV; V_fi=15; real J_fi(time) per_ms; real J_so(time) per_ms; real J_si(time) per_ms; real J_stim(time) per_ms; real p(time) dimensionless; real u_c dimensionless; u_c=0.13; real q(time) dimensionless; real u_v dimensionless; u_v=0.025; real tau_d ms; real g_fi_max mS_per_cm2; g_fi_max=8.7; real v(time) dimensionless; when(time=time.min) v=1; real tau_v_minus(time) ms; real tau_v1_minus ms; tau_v1_minus=333; real tau_v2_minus ms; tau_v2_minus=40; real tau_v_plus ms; tau_v_plus=10; real tau_0 ms; tau_0=12.5; real tau_r ms; tau_r=25; real tau_si ms; tau_si=22.22; real u_csi dimensionless; u_csi=0.85; real k dimensionless; k=10; real w(time) dimensionless; when(time=time.min) w=1; real tau_w_minus ms; tau_w_minus=65; real tau_w_plus ms; tau_w_plus=1000; real IstimStart ms; IstimStart=10; real IstimEnd ms; IstimEnd=50000; real IstimAmplitude per_ms; IstimAmplitude=-0.2; real IstimPeriod ms; IstimPeriod=1000; real IstimPulseDuration ms; IstimPulseDuration=1; // // u:time=((-1)*(J_fi+J_so+J_si+J_stim)); Vm=(V_0+u*(V_fi-V_0)); // p=(if (u q=(if (u tau_d=(Cm/g_fi_max); J_fi=((-1)*v*p*(1-u)*(u-u_c)/tau_d); // tau_v_minus=(q*tau_v1_minus+(1-q)*tau_v2_minus); v:time=((1-p)*(1-v)/tau_v_minus-p*v/tau_v_plus); // J_so=(u*(1-p)/tau_0+p/tau_r); // J_si=((-1)*w*(1+tanh(k*(u-u_csi)))/(2*tau_si)); // w:time=((1-p)*(1-w)/tau_w_minus-p*w/tau_w_plus); // J_stim=(if (((time>=IstimStart) and (time<=IstimEnd)) and ((time-IstimStart-floor((time-IstimStart)/IstimPeriod)*IstimPeriod)<=IstimPulseDuration)) IstimAmplitude else (0 per_ms)); }