|Home | About | Journals | Submit | Contact Us | Français|
The use of computational modeling to predict arrhythmia and arrhythmogensis is a relatively new field, but has nonetheless dramatically enhanced our understanding of the physiological and pathophysiological mechanisms that lead to arrhythmia. This review summarizes recent advances in the field of computational modeling approaches with a brief review of the evolution of cellular action potential models, and the incorporation of genetic mutations to understand fundamental arrhythmia mechanisms, including how simulations have revealed situation specific mechanisms leading to multiple phenotypes for the same genotype. The review then focuses on modeling drug blockade to understand how the less-than-intuitive effects some drugs have to either ameliorate or paradoxically exacerbate arrhythmia. Quantification of specific arrhythmia indicies are discussed at each spatial scale, from channel to tissue. The utility of hERG modeling to assess altered repolarization in response to drug blockade is also briefly discussed. Finally, insights gained from Ca2+ dynamical modeling and EC coupling, neurohumoral regulation of cardiac dynamics, and cell signaling pathways are also reviewed.
One of the key problems in understanding mechanisms of electrical activity in the heart, is the difficulty in predicting how perturbations that can occur at the gene, protein, cell and tissue scales alter the functioning of the whole organ. Cardiac arrhythmia is an example of altered heart function that can impair coordinated contraction that is required to maintain normal pumping and blood pressure. We now know that there are a multitude of perturbations to the heart, at every scale, that can result in cardiac arrhythmia. Examples include defects at the gene level, for example, in response to chronic drug treatment, mutations, polymorphisms, or ageing. Abnormalities such as hypertrophy can be observed in cells, at the tissue level infarct or fibrosis may occur, and deformations leading to arrhythmia may even be attributed to inherited or acquired abnormalities in organ structure. In this review, we describe some examples of simulation-based approaches to understanding cardiac dynamics in the context of arrhythmia and antiarrhythmic therapy.
Action potentials can be divided into self-oscillatory types such as pacemaker cells, and those that require an external stimulus, such as atrial and ventricular cells . Although comprehensive computational models for all the cell types have been developed and widely used, this review will focus on computational studies of the ventricle.
The ventricular myocardium displays many action potential morphologies owing to the wide variety of cell types including purkinje cells, endocardial, midmyocardial, and epicardial cells found throughout the ventricular wall. Although many action potential morphologies exist, depending on the location in the myocardium, the “classical” action potential generally has 4 phases. Phase 0 is the rapid depolarizing phase that results when Na+ channels activate and an influx of Na+ causes the membrane potential to depolarize. Phase 1 corresponds to inactivation of the Na+ channels and also outward movement of K+ ions through Ito (transient outward) currents; this contributes to the notch found in some ventricular cell types. In phase 2, a low conductance plateau phase, inward and outward ion movements are balanced by T-type and L-type Ca2+ channels and the delayed rectifier K+ channels, respectively. Phase 3 marks the final repolarization phase of the action potential, allowing the cell to return to its resting potential in phase 4. Simulations allow for important characteristics of action potential shape, morphology and duration to be quantified including the resting membrane potential (Vrest), the maximum upstroke velocity (dVm/dt), the peak overshoot value during phase 0 depolarization, and action potential duration at 50% and 90% repolarization (APD50 and APD90, respectively). These parameters are typically tracked in single cell simulations and used to compare between physiological and pathophysiological states.
It has been shown experimentally that the myocardium displays distinct cell types that are heterogeneous in ion channel expression across the ventricular wall. This heterogeneity plays a key role in the normal sequence of ventricular excitation and in arrhythmogenesis[2, 3]. Potassium currents (IKr and IKs) in particular, which are essential in ventricular repolarization, show marked differences in expression between endocardial, midmyocardial (M) and epicardial cells. By incorporating ion channel density heterogeneity into the Luo-Rudy (LR) model, Viswanathan et al.  has shown that M cells, which have reduced IKs density, display longer APD than other cell types, and these differences in multicellular models can give rise to transmural heterogeneity which may contribute to the arrhythmogenic potential in long-QT syndrome[6, 7].
Computational models have now been developed that simulate cardiac action potentials from different species including canine[8, 9], guinea pig[10-15], human[16-19], frog, and rabbit. Species specificity is important; for example, ventricular cells of higher mammalian species (canine, human, guinea pig and rabbit) have action potentials characterized by a spike and dome morphology and a well-defined plateau phase; in contrast, murine ventricular cells lack a well defined plateau phase and have a much shorter APD90[1, 23-25]. Second generation models also take into account the importance of gender, specific cell-types (including endocardial, midmyocardial, and epicardial cells) which are important for transmural heterogeneity, and include numerous pumps, exchangers, and buffers that dynamically control the intracellular ion concentrations.
These models have allowed quantification of species specific and cell-type specific functional properties of action potential morphology, response to cardiac pharmacology, and general cardiac dynamics. For example, it is known that different drugs affect ionic currents to change action potential waveforms in a species, age, and gender, and cell-type specific manner.
The use of experimental animals such as guinea pig, canine, rat and mouse have helped to build species specific models that allow comparison of simulation with experiment. However, to study phenomena such as arrhythmia, which are not necessarily the same in animal and human, requires models based on human data that recapitulate the specific characteristics of human cardiac electrodynamics. The relative paucity of human data, due both because of invasive experimental techniques, as well as insufficient methods to capture the 3D dynamics of arrhythmia initiation and propagation, makes quantitative modeling of the human heart a particularly useful tool in gaining insight into arrhythmia, mechanisms of heart disease, and physiological and pathophysiological responses to drugs and interventions.
In the past few years, models of the human cardiac ventricular myocyte have been developed from ten Tusscher[16, 18], Priebe and Beuckelmann, and Iyer that rely more on human data. The first human ventricular myocyte model published was from Priebe and Bueckelmann in 1998, and was based off of the Luo-Rudy 2 (LR2) model in which formulations for the major ionic currents were adjusted to the limited dataset available for human ventricular cells. The next substantial development in human modeling was the 2004 ten Tusscher-Noble-Noble-Panfilov model, which was based on a much larger dataset. Both models, however, still relied on human ventricular cell experiments as well as ion channel expression systems, each with their own strengths and weaknesses.
Human models (and others) have increased greatly in complexity and require sufficient computational power for higher dimension simulations. Both the ten Tusscher and Priebe and Beuckelmann models are based off of the Hodgkin-Huxley (HH) formalism and have approximately 16 variables, while the Iyer model is Markov-based and has approximately 67 variables. These differences are important to consider with regard to application; large-scale spatial simulations might require the use of less complex models such as the ten Tusscher or Priebe and Bueckelmann models that are computationally more efficient. Conversely, ion channel pharmacokinetic studies requiring the resolution of simulating single ion channel states might require the Markov formalism found in the Iyer model. A detailed comparison of models for human ventricular cells and tissues can be found in .
An important step in understanding how arrhythmias are initiated, sustained and terminated is to develop quantitative indices of “arrhythmia vulnerability”, to identify the underlying “vulnerable parameter” and to ultimately reveal the ionic mechanism of the vulnerable parameter. Examples of vulnerable parameters of an arrhythmia include effective refractory period, reduced excitability, dispersion of repolarization and conduction block. For example, refractory cells and tissue are inexcitable, due to insufficient time for recovery of Na+ channels from time- and voltage-dependent inactivation, a condition that can be exacerbated by mutations or Na+ channel block by drugs. Dispersion of repolarization can occur when premature impulses propagate through cardiac tissue.
The dynamical factors that govern such spatial gradients are APD restitution (the shortening of APD as a function of shortening diastolic interval (DI)) and conduction velocity (CV) restitution (the slowing of CV as a function of shortening DI) – these two factors cause wave dynamics to change as a function of space for premature impulses and set the stage for functional block and arrhythmia initiation. These factors are also quantitative indices for arrhythmia vulnerability.
The primary tissue-level cause of arrhythmia initiation is conduction block. Conduction block occurs when the spatial propagation of an action potential terminates prematurely. Conduction block sets the stage for reentry, in which the action potential propagates around the site of block and later activates it in the reverse direction, thereby starting a reentrant circuit.
Because arrhythmia is fundamentally a spatiotemporal phenomena, markers of arrhythmia vulnerability change at different spatial scales (channel, cell, tissue), and some properties that seem to be antiarrhythmic at one level (for example increased refractory period in single cells), are in fact proarrhythmic at the tissue level (see section on vulnerable window). Mathematical modeling and simulation can therefore track and quantify markers of vulnerability at each spatial scale and determine their relative importance in initiating arrhythmia. We begin with markers at the cellular level.
The Luo-Rudy (LR) models, a paradigm of cardiac ventricular action potential modeling, were among the first to study mechanisms of arrhythmogenesis in simulated ventricular myocytes. In 1991, the first Luo-Rudy model (LR1) newly formulated 4 ionic currents: INa, the fast Na+ current, which included fast and slow inactivation gating; IK, the delayed rectifier current; IK1, the inward rectifier K+ current; and a plateau current, IKp . The reformulation of INa with a much larger channel conductance, allowed the model to reproduce action potential upstroke velocity typically seen in single cardiac myocytes (dVm/dtmax≈ 400 V/sec). LR1 was most significant for reproducing cardiac dynamics that were dependent on [K+]o including APD and resting membrane potential. By incorporating [K+]o as a dynamic variable, [K+]o – dependent phenomena such as responses to premature stimuli and periodic stimulation reproduced monotonic Wenckeback patterns and alternans at normal [K+]o, and nonmonotonic Wenckeback periodicities, aperiodic patterns, and enhanced supernormal excitability at low [K+]o . These observations were consistent with experimental observations[33, 34], and suggested that intrinsic properties of single cardiac cells such as Wenckebach periodicities and aperiodic responses to periodic stimulation might be important in arrhythmogenesis.
In the 1994 Luo-Rudy papers, a description of cell processes that regulated intracellular ion concentrations of Na+, K+ and Ca2+ were included. Ca2+ dependent processes were described in detail and new formulations for an L-type Ca2+ channel, a Na+/Ca2+ exchanger, Na+/K+ pump, two compartment sarcoplasmic reticulum (SR) Ca2+ uptake and release, and a nonspecific Ca2+ activated membrane current were added. Release of Ca2+ from the SR could be triggered by both Ca2+– induced – Ca2+ release and spontaneous release. Buffering of Ca2+ was implemented with both troponin and calmodulin for the cytoplasm and calsequestrin for the SR.
The inclusion of Ca2+ dynamics, in conjunction with dynamic calculations of [Na+]i and [K+]i allowed for simulations of triggered activity, EADs, and DADs in the accompanying article. Two different mechanisms of EADs were proposed that depended on the phase of the action potential in which they occurred (phase 2, “plateau”; and phase 3 EADs). For the first time, modeling studies were able to quantitatively define mechanisms of afterdepolarizations and triggered activity, which were the result of a complex interplay between many interacting processes with regard to Ca2+ handling.
The Luo-Rudy models have also been used as tools to study mutation induced arrhythmia. For example, by incorporating a model of a mutant Na+ channel harboring long-QT3 (LQT3), simulation studies showed that the mutant phenotype of LQT3 (aberrant Na+ channel gating leading to a small, sustained inward current), is not directly responsible for triggering the phase 2 EADs commonly seen in experiments. Rather, the sustained current sets up favorable conditions (prolonged membrane depolarization at plateau potentials) that allow sufficient time for L-type Ca2+ channels (LTCCs) to recover and reactivate, providing the actual trigger for EADs.
Ventricular fibrillation is thought to arise from a dynamic instability in cardiac tissue that can cause APD alternans (an instability in cellular dynamics that results in a beat-to-beat alternation in APD – short, long, short, long etc.) leading to fragmentation of spiral waves and fibrillation-like excitation patterns. Mathematically, APD restitution describes the relationship between the APD and the preceding diastolic interval (DI):
where CL is the cycle length of the static pacing regime.
The reduction of APD as a function of shortening diastolic interval (DI) is the APD restitution relation and forms the basis for the restitution hypothesis which states that a steep restitution curve with a slope >1 is the criterion for the onset of alternans, APD instability and ventricular fibrillation. Conversely, flattening the APD restitution curve inhibits alternans and subsequent conduction blocks, thereby preventing fibrillation. The graphical method of Nolasco and Dalen mathematically shows that when the slope of the APD restitution curve is >1, the APD fails to converge with periodic beating and gives way to alternans.
When restitution is steep, spatial dispersion of refractoriness can occur for single premature beats or for a series of beats. Alternans can also occur spatially, where cells in different regions alternate in APD. During alternans, dispersion of refractoriness can be monotonic (i.e., consistently increasing or decreasing) or discontinuous (i.e., a phase reversal of the gradient that occurs at a point known as a “node”). An example of a discontinuous case is discordant alternans, in which a cardiac wave has a short APD in one (or more) region(s) of tissue but a long APD in another region.
Many experimental studies have measured restitution curves which show not only a broad range of restitution slopes[38-40], but also spatial variability within the myocardium. APD restitution in the human heart has been measured both close to 1[41, 42], and considerably steeper – from 1.28 in the apex to 3.79 in the outflow track, with a positive correlation between arrhythmia inducibility, steepness and dispersion of repolarization.
While experiments[44, 45] and simulations[46, 47] have both shown that ventricular fibrillation can be suppressed by flattening the restitution curve pharmacologically and preventing wave-break instabilities, other experiments have shown preparations with APD restitution curve slopes >1, which have exhibited stable behavior without alternans development, even at short cycle lengths where slopes are steepest[40, 48]. This suggests that steep restitution may not be the singular criterion for APD alternans.
In fact, while modeling studies have shown that steep restitution is necessary for alternans, APD restitution slope >1 may be an oversimplification that only holds for very simple models. Other factors such as short-term cardiac memory, electrontonic interactions between cells, conduction velocity (CV) restitution, the range of diastolic intervals (DIs) over which restitution is steeper than 1, and the range of DIs visited during spiral wave rotation are also important determinants of alternans and spiral wave breakup[16, 49-54]. Notably, the study by Cherry and Fenton use modeling to explain under what conditions electrotonic effects and memory can suppress alternans and conduction block despite an APD restitution curve slope much greater than 1.
Lastly, studies suggest that conduction velocity (CV) restitution (the dependence of the CV on preceding DI), together with APD restitution, also play an important role in electrical dynamics of cardiac tissue[45, 55-57]. Incorporating multiple electrical dynamical mechanisms to explain propensity for arrhythmia have revealed that simulations and modeling studies are fundamental to probing how mutations and drugs alter tissue dynamics and lead to the generation of arrhythmia.
One of the challenges of understanding arrhythmia and designing effective antiarrhythmic drugs is the linking of antiarrhythmic action from electrophysiological experiments to tissue and organ level effects. For a long time, it was unclear why the CAST trials[58, 59], one the largest placebo controlled trials of antiarrhythmic drugs, failed so profoundly. Drugs that were shown to be antiarrhythmic at the single cell level paradoxically increased mortality by 2-3x as compared to placebo[58, 59]. As stated by Sanderson, “In few specialties of medicine are new promising drugs shown to be so much inferior to placebo, and even worse, to increase mortality”.
Briefly, potentially fatal cardiac arrhythmias can arise from unexpected stimuli or premature ventricular contractions (PVCs). If the PVC is sufficiently strong to excite and recruit neighboring cells, a continuous wave or wave fragment will form, depending on the excitability of those neighbors. A continuous wave that propagates in both directions will eventually collide with itself and extinguish; however, a wave fragment that propagates in some directions, but fails in others can lead to wavelets and spiral waves. This is one purported mechanism of reentry, and self-sustaining oscillations producing a functional reentrant arrhythmia can lead to loss of heart pump function and death.
For a long time, the prevailing theory in treating cardiac arrhythmia was drugs that increased refractoriness to PVCs (or ectopic beats), would reduce the incidence of arrhythmias leading to sudden cardiac death. This can be accomplished by reducing excitability by Na+ channel blockade and prolonging the refractory period, or by K+ channel blockade, which acts to prolong the action potential duration. Interestingly, although prolongation of the recovery process does result in PVC suppression, it simultaneously increases the duration of the vulnerable window (Figure 1) – a time period where a PVC can trigger a reentrant arrhythmia – (a proarrhythmic effect), exposing unsuppressed PVCs to more vulnerable tissue. So, while ion channel blockade decreases the number of PVCs, it increases the likelihood that unsupressed PVCs will initiate an arrhythmia.
Ferris et al. and Wiggers and Wegria observed that ventricular fibrillation could be readily induced by stimulating the heart in an interval they called the vulnerable phase. Later, it was shown that stimulation within the vulnerable region would lead to spiral waves and ventricular fibrillation.
Computational modeling is an extraordinarily useful tool to define the parameter space of vulnerability in a simple one-dimensional (1D) system that can then be extended and applied to spiral wave induction in twoand three-dimensions. Consider a homogeneous fiber stimulated (S1) at cell 0 for various pacing cycle lengths (PCLs) until steady state is reached. Application of a premature (S2) stimulus is then applied in the middle of the fiber at varying inter-stimulus intervals (relative to S1). If the interval is too short, the stimulus hits refractory tissue from the preceding wave and fails to propagate in either direction (bidirectional block). Continuing to increase the inter-stimulus interval will eventually initiate a wave that propagates retrogradely, but not in the antegrade direction; this marks the transition from the absolute refractory period to the “most premature beat (MPB)” that causes retrograde conduction. A further increase in the inter-stimulus interval will cause a second transition, or “least premature beat (LPB)”. This will delineate the transition from unidirectional retrograde conduction to bidirectional conduction (the presence of both an antegradely- and retrogradely- propagating wave). This time interval (tLPB – tMPB) is denoted the vulnerable window (VW) for retrograde conduction (1D) and analogously spiral wave initiation in 2D. Further details can be found in Starmer et al.[61, 62, 65, 66].
Vulnerability analysis clearly elucidates the often less-than-intuitive link between functional levels of arrhythmia suppression and helps to answer the paradox put forth earlier. The same mechanism (increased refractoriness) that is antiarrhythmic in single cells is in fact proarrhythmic in higher dimensions! For a more in-depth review of the cardiac vulnerable window, see[61, 62, 65-67].
Numerous mutations in cardiac ion channels have been causally linked to a number of diverse clinical syndromes including Long-QT Syndrome (LQTs), Brugada Syndrome (BrS), and Isolated Cardiac Conduction Disorder (ICCD) (Figure 2). Because the literature on the hundreds of identified mutations and dozens of computational studies investigating their mechanisms is so vast, we will focus this review on a survey of computational studies investigating mechanisms of Na+ channel linked inherited arrhythmias.
Mutations in cardiac Na+ channels that underlie LQTS tend to disrupt the critical property of channel inactivation. The phenotypic manifestation of LQTS is a prolongation of the QT interval on the ECG. In general, Na+ channel linked LQTS stems from mutation-induced disruption of channel inactivation, as was originally identified in the ΔKPQ mutation, a three amino acid deletion in the intracellular linker between domains III and IV of NaV1.5. This motif is known to be critical for fast inactivation of the channel and the mutation results in persistent non-inactivating current. The non-inactivating component of INa, acts to prolong the plateau of the action potential and may allow for the development of arrhythmogenic triggered activity, referred to as early afterdepolarizations (EADs)[70-72].
While ΔKPQ is one example of altered gating, evidence suggests that mutation induced gain of function in cardiac INa can exist in at least three distinct forms. Computational approaches have been implemented to show that any of these altered forms of gating can lead to the observed long-QT phenotype. Channel bursting, characterized by long channel open times is apparently due to transient inactivation failure - as in ΔKPQ - that underlies sustained Na+ channel activity over the plateau voltage range. Computational investigation has shown the extent of bursting that is required to cause triggered activity and arrhythmia[70, 73]. A second abnormality in Na+ channel gating that has been computational tested and linked to phenotype results from steady state channel reopening called window current[68, 74], because reopening occurs over voltage ranges for which steady-state inactivation and activation overlap. A third original mechanism was demonstrated in channels containing the I1768V mutation, which does not result in an obvious gain of channel function. However, simulations, in conjunction with experiments revealed that under non-equilibrium conditions during repolarization, channel reopening results from faster recovery from inactivation at membrane potentials that facilitate the activation transition. Mutation induced faster recovery from inactivation results in channels that reopen during repolarization and the resulting current amplitude rivals that of bursting channels. Simulations have demonstrated that late current due to channels reopening causes severe prolongation of the AP plateau and arrhythmic triggers.
In contrast to LQTS associated mutations that cause a gain of Na+ channel function, mutations that have been identified in patients with congenital forms of Burgada syndrome (BrS) result in a loss of Na+ channel function, such as the Y1795H mutation in the cardiac C-terminus[75-77]. BrS is an arrhythmic syndrome characterized by right bundle branch block, ST segment elevation on the ECG and ventricular tachyarrhythmias that result in sudden death. The reduction in INa has been shown to occur by several mechanisms in BrS, including reduced rates of recovery from inactivation, faster inactivation subsequent to channel opening and protein trafficking defects[79-82].
One proposed mechanism underlying the ECG manifestations in BrS is that heterogeneity of repolarization exists across the right ventricular myocardial wall. Repolarization in the epicardial cell layer is thought to be more sensitive to suppression of the peak Na+ current due to a more prominent transient outward K+ current (Ito). Mutations in SCN5A that reduce the peak INa could selectively hasten the repolarization in the epicardium, leading to a loss of the epicardial action potential plateau phase and an increased propensity towards reentrant arrhythmias. This theory may explain the ST-segment elevation resulting from a dispersion of action potential plateau potentials observed in patients with BrS. Several computationally based investigations have suggested the plausibility of this mechanism[84-86].
Another proposed mechanism underlying the Brugada syndrome phenotype is that loss of Na+ current, observed as a depolarizing shift of the channel activation curve, may result in sufficiently slow conduction that some areas of the myocardium are undergoing repolarization while others are still depolarizing. Indeed, a simulation- based investigation of this mechanism was investigated and shown to be plausible and account for the observed ST-segment elevation in these patients[87-89].
The model simulations by Tan et al.  investigated the isolated cardiac conduction disease mutation G514C. The showed in their simulations that a marked decrease in conduction velocity occurred as activation gating is right shifted. There was also a concomitant increase in the required strength of stimulus current to excite coupled tissue. The investigation revealed successful reversal of the pathological phenotype upon glucocorticoid steroid treatment and restoration of electrical propagation.
The relationship between genetic mutations and clinical syndromes is complex, as the revelation of novel mutations suggests paradoxical phenotypic overlap or exclusivity. Multiple loci in the cardiac Na+ channel have been identified where the same mutation can result in different disease phenotypes. An insertion of an aspartic acid residue (1795insD) in the carboxy terminus of NaV1.5 can result in either BrS or LQTS. Additionally, mutation of the same residue to a histidine (Y1795H) or cysteine (Y1795C) results in BrS and LQTS, respectively, indicating the proximal region of the C-terminus as a potentially important structure in sodium channel function.
The deletion of a lysine (*K1500) in the III-IV linker of NaV1.5 is associated with BrS, LQTs and ICCD[79, 80]. LQTS is typically associated with gain-of-function Na+ channel mutations while BrS and ICCD are typically associated with loss-of-function, resulting in reduced INa. How can a single mutation simultaneously result in seemingly paradoxical syndromes (i.e. gain of function LQTS and loss of function BrS)?
Simulation studies have suggested that one explanation may stem from the intrinsic heterogeneity of the underlying myocardial substrate with which mutant Na+ channels interact. The ventricular myocardium is comprised of at least three distinct cell types (epicardial, midmyocardial (M) and endocardial cells), which exhibit distinct electrophysiological properties. Epicardial cells display a characteristic spike and dome morphology due to large transient outward K+ current (Ito), and short action potential duration (APD) resulting from a high density of the slowly activating component of the delayed rectifier K+ current IKs . Mutations that act to reduce INa in the presence of large repolarizing currents (Ito and IKs), may result in premature plateau repolarization (BrS phenotype) and APs with distinctive triangular morphology or in action potentials with prominent coved domes[83, 84, 86, 93, 94]. Ito and IKs are smaller in M-cells, and are unable to overwhelm the mutation induced reduced INa. Selective loss of the AP plateau in epicaridal cells, results in dispersion of plateau potentials across the ventricular wall. This gradient generates ST segment elevation on the ECG, which is a diagnostic indicator of BrS[83, 95]. Clinically, ST segment elevation is observed in right precordial leads of BrS patients, consistent with the large Ito density in right ventricular epicardium[83, 86]. In M-cells the non-inactivating component of INa, in the presence of smaller repolarizing currents, acts to prolong the plateau of the action potential, and may allow for the development of arrhythmogenic early after depolarizations (LQT phenotype). APD prolongation is reflected in a prolonged QT interval on the ECG, indicative of the LQT Syndrome.
There are multiple recent examples of the importance of pharmacogenetic considerations in drug therapy. For example, the recent discovery of a common polymorphism present in 13% of African Americans that predisposes carriers to drug induced arrhythmia. Simulations revealed the mechanism of the increased susceptibility of carriers to hERG block. Another study examined the effects of the Na+ channel open state blocker mexiletine and the inactivation state blocker lidocaine on LQT3 associated ΔKPQ mutant channels. The simulations suggested that open state blockers preferentially reduce the arrhythmogenic persistent current that is the hallmark of LQT3 mutations in cardiac Na+ channels. Moreover, this preferential block is achieved at relatively low doses of the drug that have little effect on tissue conduction. The study suggests that low doses of mexiletine will effectively normalize action potential duration at slow heart rates in the presence of the LQT associated ΔKPQ mutation, while application of lidocaine will not. Moreover, the model suggests that the effective dose of mexiletine required to reverse the LQT phenotype has minimal effects on peak current, and therefore on the action potential upstroke velocity, a primary determinant of conduction velocity. Lidocaine, on the other hand, causes large reduction of peak current. The results suggest Na+ channel open state block as a genotype-driven potential therapeutic target.
Computational modeling can also be used to assess off target effects of drug toxicity due to hERG inhibition. It has been suggested that approximately 40% of all new drugs interfere with the repolarization process, which can lead to malignant arrhythmia and cause a drug to be pulled from the market. Accurate prediction of drug / hERG interaction would thus be beneficial for the drug development cycle, by discarding candidate compounds much earlier in the development process. However, current markers to assess drug-induced arrhythmia, QT interval and IKr inhibition, both have limitations in predicting clinical outcome. Computation modeling is powerful tool and may be one solution to this problem, as it can follow makers of arrhythmia from channel level and specific ionic currents, to clinical level (e.g. QT interval). In a recent study, computational modeling was used to understand mechanism(s) responsible for a hERG channel opener's (NS1643) antiarrhythmic effect both in normal, and hypokalemic conditions. Their results interestingly conclude that an increased hERG conductance and a depolarizing shift of the inactivation curve both contribute to the antiarrhythmic actions of NS1643, with relative effects dependent on external K+ concentration.
In another recent study of heteromeric hERG 1α/1β channels underlying cardiac IKr, computational kinetic modeling was used to determine whether the faster apparent rates of activation and recovery from inactivation could quantitatively account for the increased current amplitude of hERG 1α/1β channels and to explain a reduction in efficacy of drug block in the heteromeric hERG channels. A model investigation revealed that an alteration in kinetics is sufficient to account for reduced rectification and a concomitant increase in current amplitude during depolarization. The simulations also showed reduced E-4031 sensitivity for the hERG 1α/1β heteromer, due to the absence of additional E-4031 blocked states in an additional mode of gating unique to the homomeric channel. Finally, this study showed that for increasing doses of E-4031, APD90 becomes longer for both hERG 1α/1β and hERG 1α. Action potential simulations revealed that pause-induced early afterdepolarizations are evoked selectively for hERG 1α. Taken together with the fact that excessive APD prolongation can cause arrhythmia, these results suggest that hERG blocking drugs are more arrhythmogenic for hERG 1α than they are for hERG 1α/1β.
Accurate modeling of calcium dynamics is one of the most complex areas of cardiac action potential simulation. It is intrinsically complicated and highly integrated with multiple components including the L- and T-type Ca2+ current, ryanodine receptors, Na+/Ca2+ exchangers and SERCA pumps, Ca2+ binding proteins and buffering with calmodulin and troponin, and multiple sarcoplasmic reticulum (SR) subspaces. Detailed formulations of Ca2+ dynamics serve as the link between electrical depolarization of the cellular membrane and activation of the contractile apparatus, so-called “EC coupling”. Part of the difficulty in building accurate models of the cardiac Ca2+ handling system stems from the seemingly paradoxical requirement of a high gain and graded response system that is under tight control of ICa, as well as multiple feedback and feedforward loops [2, 22].
EC coupling has a general mechanism as follows: Ca2+ influx into the sarcolemma induces a secondary Ca2+ release from the sarcoplasmic reticulum (so-called CICR (Ca2+ – induced – Ca2+ release)) via ryanodine receptors. Ryanodine receptors sit in a restricted didactic space and sense local Ca2+ concentration. Ca2+ influx from the L-type Ca2+ channel (LTCC) and release from the SR form a negative feedback loop to induce inactivation of the LTCC and reduce further influx. The released Ca2+ is then buffered by calmodulin and troponin, which generates force via myofilament activation. Most Ca2+ is resequestered via the SERCA pump each heartbeat with the remainder extruded via the Na+/Ca2+ exchanger[100, 101].
In the early 1990s, it was realized that to fully capture the cyclical nature and essential features of EC coupling, modeling of Ca2+ dynamics needed to incorporate subcellular compartments. Nordin et al. developed a 3 compartment model that relied on gradient diffusion laws and empirical adjustment. More complex models have since been developed that incorporate subsarcolemmal microdomains for local Ca2+ sensing with three levels of organization: the cleft for LTCC and ryanodine receptor interactions; the subsarcolemmal space for interactions between transporters and their corresponding ions; and the cytosolic space where force is generated by interaction between Ca2+ and myofilaments.
These models have provided insight into EC coupling that is far more complex and integrated than intuition alone could predict. For example, the requirement of minute volume subspaces for Ca2+-induced-Ca2+release, and the interaction between ryanodine receptors and the LTCC, which includes the effects of surface charge and the electric field in the cardiac diad. The Soeller and Cannel model study showed that the surface charge has significant effects on the time course of calcium movements in the diad and causes local Ca2+ changes to occur more rapidly than predicted from models that consider sarcolemmal Ca2+ binding without field effects.
However, as a result of Ca2+ complexity, arrhythmias resulting from aberrant Ca2+ dynamics form a large subset of congenital and acquired arrhythmias. Most often, these arrhythmias are the result of Ca2+ overload, which can lead to delayed afterdepolarizations (DADs), premature ectopic beats (extrasystoles), monomorphic and polymorphic ventricular tachycardia, bidirectional tachycardia, electrical alternans, atrial fibrillation, and AV nodal arrhythmias. Both from a physiological and pathophysiological standpoint, accurate models of Ca2+ handling dynamics and EC coupling have been important in integrating multiple interacting processes to suggest mechanisms leading to arrhythmia initiation. For example, it was modeling that first implicated the LTCC in the genesis of EADs secondary to persistent Na+ due to mutations.
Another example is the modeling of spontaneous Ca2+ release during SR overload, a feature common to all mammalian cardiac tissue. Oscillatory Ca2+ release and ectopic beating can be induced through inhibition of the Na+/K+-ATPase, which produces Na+ overload and subsequent Ca2+ overload. Nordin was able to simulate the production of DADs and ectopic beating through a similar mechanism in a ventricular cell model which relied on high gain in the SR-release mechanism as SR load increased during the overload. The LR2 model also demonstrated ectopic beating, spontaneous release and DADs, but implemented a heuristic mechanism that triggers Ca2+ release after a set fraction of Ca2+ bound to calsequestrin exceeded a threshold[12, 100].
Because it is only recently that mutations directly linked to specific defects in myocardial Ca2+ handling have been elucidated, the field of inherited Ca2+ dependent arrhythmia syndromes, such as mutations in the ryanodine receptor (RyR2) leading to catecholaminergic polymorphic ventricular tachycardia (CPVT)  and arrhythmogenic right ventricular cardiomyopathy type 2, will be a fruitful avenue of research for mathematical modeling, such as has been the case with inherited Na+ channelopathies. In a recent study, Iyer et al. simulated mutant RyR2 in a model of CPVT; the model proposes that hyperactive, “leaky” receptors characteristic of reduced FKBP12.6 function may be centrally involved in triggering DADs, and may provide a plausible link between defective RyR2 gating and CPVT. Ca2+ handling and signaling leading to arrhythmia is decidedly a complex and in-depth topic, the reader is referred to more detailed reviews of Ca2+ handling in cardiac cells, and myocardial Ca2+ signaling and arrhythmia pathogenesis.
Cell signaling networks are indispensible for many cellular functions and form a tightly interwoven network that coordinates and controls aspects of gene regulation, metabolism, chemotaxis, electrical activity and contractile behavior. Because of the multiple levels of interdependencies, feedback and feedforward loops, a quantitative understanding of cell signaling both in relation to physiological and pathophysiological conditions remains a fundamental challenge to both experimentalists and modelers.
Most models of the cardiac ventricular myocyte assume steady state levels of neurohumoral response of adrenergic and other cell stimuli, and thus neglect cell-signaling pathways in action potential simulations. However, it is clear that cell signaling has profound effects on many cellular process that regulate cardiac electrodynamics; for example β-adrenergic modulation of IKS, indirect modulation of Ca2+ dynamics, and the L-type Ca2+ channel, to name a few.
Briefly, the β-adrenergic signaling cascade relies on catecholamines to coordinate control of contractility, metabolism and gene regulation[111, 112]. β receptors couple with the GPCR subunit GS to stimulate adneylyl cyclase, which synthesizes cAMP to act as a second messenger, promoting dissociation of PKA. PKA then phosphorylates and alters function of multiple target proteins and processes including the L-type Ca2+ channel, phospholamban, ryanodine receptors, INa, IKs, Na+/K+-ATPase, sarcolemmal Ca2+-ATPase, Troponin I, crossbridge cycling rates, and glycolysis[22, 113].
A problem arises with the assumption of either steady state levels of neurohumoral response or bulk cytoplasmic concentrations of second messengers in signaling cascades. The need for models to incorporate microdomain-level processes is evidenced by the Iancu model of cAMP signaling in cardiac myocytes. This model tested whether the effects of β1-adrenergic receptor and M2 muscarinic activation involve compartmentation of cAMP. While the bulk cAMP concentration inside the myocyte is approximately 1μM[114-116], the binding affinity for PKA is on the order of 100 – 300nM[114, 117, 118]. By segregating cAMP within the cell into two submembrane microdomains (caveolar and extra-caveolar), and a bulk cytoplasmic domain, the model was able to show that while bulk basal cAMP levels are too high to appropriately regulate PKA, caveolar cAMP varies in a more appropriate range for regulation (~100nM - ~2μM). The same group later validated many of the hypotheses generated by the computational model in a subsequent experimental study.
With regard to disease processes leading to arrhythmia, it has been recently shown that the β-adrenergic signaling cascade involves the association of PKA and PP1 with the C-terminal domain of KCNQ1 (the α-subunit of IKs). β-adrenergic modulation of IK1 in LQT1 patients upsets the delicate balance of repolarization and can lead to arrhythmia; this is the rationale for β-blockade pharmacotherapy for this subset of LQT patients. Interestingly, simulation studies have questioned the use of such therapy as a general treatment strategy for all LQT patients, and may actually be proarrhythmic in LQT3 patients. Other diseases, such as catecholaminergic polymorphic ventricular tachycardia, are linked to abnormal Ca2+ cycling due to increased β-adrenergic tone[122, 123].
Saucerman et al. has used computational approaches to link β-adrenergic signaling with excitation contraction coupling and genetic mutations[108, 111, 124] (see Figure 3). In one simulation study, Saucerman et al. found that a known gene mutation (KCNQ1-G589D) associated with LQT syndrome did not alone prolong the QT interval; rather it was only when the model was coupled to β-adrenergic stimulation, did the mutation induce QT prolongation and transient afterdepolarizations. When implemented into a 3D rabbit ventricular wedge model, the mutant further elevated transmural dispersion of repolarization and created T-wave abnormalities on simulated ECGs, all findings suggestive of a mutation-induced disruption of β-adrenergic stimulation and electrophysiology leading to known triggers of cardiac arrhythmia.
This model and others[21, 108, 124] highlight that the link between neurohumoral regulation, Ca2+ cycling and cardiac electrodynamics is highly complex, with multiple levels of integration and subcellular localization. Most of these second messengers are not free to diffuse throughout the cell, but rather form compartmentalized networks, which allow for differential responses to multiple stimuli, control cross-talk, and coordinate distinct phosphorylation targets[22, 125-127].
In sum, prediction of cellular responses to any aberration in cell-signaling cascades or neurohumoral regulation by disease, drugs and other interventions would be virtually impossible without the use of predictive modeling. “Third generation” cardiac ventricular myocyte models will be notable for inclusion of complex signaling pathways and will most likely rely on a systems level prospective to understand how perturbations in any of these networks and components will affect cardiac dynamics.
Simulation of cardiac action potentials is important for understanding the effects of perturbations at the ion channel scale such as mutant channel aberrant gating and drug / channel interaction, but is also important as the foundation for higher dimensional simulations that integrate knowledge at different spatial scales into a framework for understanding physiology and pathophysiology of cardiac arrhythmia and heart disease.
The field of cardiac simulation and experimental electrophysiology have worked together for the past 50 years to advance our understanding of cardiac electrophysiology through iterative simulation and experiment, and vice versa. Noble and Rudy note that as opposed to other fields in science and engineering, such as the aerospace and automotive industry, modeling of complex biological phenomena has been slow to take root. However, simulation has the potential to compliment and validate experiments and allow for prediction before suitable experimental technology becomes available to validate these results.
Cardiac action potential modeling has been successful as an approach to understanding cardiac arrhythmia and other physiological mechanisms because they are able to integrate multiple functional levels (gene, cell, tissue, organ, electrocardiogram) and the complex interplay between them, that allow for arrhythmia initiation. Traditional single-scale approaches typically fail to reveal how disruption in protein function due to mutations and/or drugs and, consequently, through complex interactions and behaviors of cells, lead to loss of synchronization and arrhythmogenic rhythms in tissue to cause failure of coordinated contraction. Integration of different spatial scales allows the computational biologist to follow key parameters over multiple spatial scales and identify interacting components.
Sufficiently detailed models of the cardiac action potential now exist at the single cell level and have allowed for tremendous insight into ion channel biophysics, cell-cell interactions, action potential propagation and other cardiac dynamics. However, these models are complex and require the combined efforts of physicians, computational biologists, biophysicists and computer scientists to integrate knowledge across disciplines. The challenge ahead remains advancing computational power to be able to integrate sufficiently detailed models into 2D tissues, whole organs or even virtual torsos to study these phenomena. Integrated understanding from genes to electrocardiogram is slowly becoming a reality .
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.