Cardiac arrhythmias represent a leading cause of morbidity and mortality in developed countries. In spite of intense research, the mechanisms of generation, maintenance and termination of theses arrhythmias are still not clearly understood. The ability to predict, prevent and treat these arrhythmias remains a major scientific challenge. In the last years, new technologies have produced a huge amount of information at different levels: subcellular, cellular, tissue, organ and system. However, the knowledge of each component is not sufficient by itself to understand its behavior and an integrative approach, which should take into account the complex interactions between all components is needed. This understanding is necessary to improve prevention and treatment of cardiac arrhythmias.
Traditionally, cardiac research has been based on experimental and clinical studies. However, in the last decades, a quantitative understanding of the biophysical and biochemical properties of individual cardiac myocytes and of the anatomical structure of the heart has permitted the development of computational models of single myocytes, portions of tissue and even the whole heart. These models have opened a new approach in studying the complex mechanisms underlying cardiac arrhythmias by means of computer simulations.
Integrated and multiscale models of the heart are improving very fast, as a large body of research is being made in this field. These models integrate, with a great degree of detail, the different levels of the organ, from the genetic characteristics of membrane ionic channels, to the electrophysiological characteristics of cardiac cells and to the anatomical structure of the different cardiac tissues. In this chapter, the basis of cardiac multiscale modeling will be reviewed and recent and specific case studies for the different levels will be described focusing on the understanding of cardiac arrhythmias during myocardial ischemia and drug treatment.
2. Modeling cardiac ion channels activity
Comprehensive models of the electrical activity of cardiac ion channels were formulated for the first time in 1952 (Hodgkin & Huxley, 1952). These models were based on the Hodgkin-Huxley mathematical formalism, which allows the reproduction of macroscopic ionic currents. This approach computes the conductance of a certain ion channel as a function of the open probability of each gate of the channel and the maximum conductance of the membrane for that particular ion. In additon, gate transitions from closed to open and vice versa follow a first order voltage dependent behavior at a rate that is independent of the remaining gate states.
The improvements in the experimental techniques employed to measure and characterize the ionic currents responsible for the cardiac action potential have provoked a paralell evolution in the ion channel models, which have acquired an increasing complexity and realism. Two examples are the separation of certain currents into its components, like the slow and rapid delayed rectifier potassimum currents, and the discovery of ion current dependencies on electrolyte concentrations. Importantly, experimental data have evidenced the existence of state dependent transitions of ion channels that Hodgkin-Huxley formalism fails to reproduce. This fact has favored the development of Markov models for ion channels, which are based on the assumption that the transitions between channel states depend on the present conformation of the channel and not on its past history. Furthermore, recent developments in molecular biology and in the genetics of ion channels have made possible the characterization of the electrophysiological properties of a certain mutated ion channel (Tomaselli et al., 1995), allowing the development of its mathematical model using Markov models. Although Markov models can provide a more realistic formulation of ion channels, Hodgkin-Huxley models are also currently used because their computational requirements are significantly less demanding. Both formalisms are used to model ion channel-drug interactions and pathological situations.
2.1. Modeling myocardial ischemia at the ion channel level
The ionic current and action potential models developed in the last decades make it possible to simulate the electrical activity of cardiac cells both in normal and abnormal conditions. Among the latter, acute myocardial ischemia is one of the most threatening abnormalities that can occur in cardiac tissue because it can easily trigger malignant and potentially mortal reentrant tachyarrhythmias such as ventricular tachycardia (VT) and ventricular fibrillation (VF) (Janse & Wit, 1989). Indeed, VF subsequent to ischemia is considered as the first cause of mortality in Europe, the USA and a significant part of Asia (Ross, 1999).
Myocardial ischemia occurs when a coronary artery is occluded. When this happens, a certain zone of the myocardium (formed by the cells that depended on the occluded artery for their blood and nutrient supply) begins to suffer metabolic and electrophysiological changes. Specifically, K+ ions begin to rapidly accumulate in the extracellular medium (a phenomemon termed hyperkalemia), the media (both intra and extracellular) become acidotic (with a drop in pH from 7.2-7.4 to 6.2-6.4 in 10-15 minutes), and oxygen levels decline (hypoxia). Hypoxia provokes in turn a depletion in intracellular ATP levels and a concomitant increase in ADP levels, which results in the partial activation of a specific type of K+ current called ATP-sensitive potassium current (IK(ATP)) (Noma, 1983).
To properly and predictively simulate the cardiac arrhythmias arising from myocardial ischemia, the model must be of a multi-scale nature and ischemia has to be modeled in the ion channel level, the cell level and the tissue level. While tissue modeling of ischemia will be adressed further in this chapter, ion channel ischemia modeling should be done as follows. First, hypoxia must be included in the model through its effects in IK(ATP) current activation (see below). Second, hyperkalemia must be simulated by increasing [K+]o levels (Rodriguez et al., 2002). Elevated [K+]o in turn alters the conductance of a number of sarcolemmal channels, such as IK1 or IKr, and these effects are included in most ion channel models published to date. Third, acidosis must be modeled through its effects in the maximum conductances and kinetic properties of the Na+ and Ca2+ inward currents (Ferrero, Jr. et al., 2003; Shaw & Rudy, 1997a).
The IK(ATP) current, which is normally dormant in normoxic situations but is partially (though only slightly) activated in acute myocardial ischemia, is carried by K+ ions (Noma, 1983), shows voltage-dependent inward rectification (but is otherwise non-voltage dependent) caused by Mg2+ and Na+ ions (Horie et al., 1987), and is gated by intracellular ATP with intracellular ADP playing a modulating role in the ATP-mediated deactivation (Weiss et al., 1992). All these facts have to be present in an ion channel model of the current. In 1995, Ferrero et al. published the first comprehensive model for the myocardial IK(ATP) current (Ferrero, Jr. et al., 1996), which has been further improved more recently (Michailova et al., 2005).
The mathematical model for the IK(ATP) current includes all the aforementioned effects, which are represented in the following equation:
where σ is the channel density, g0 the unitary conductance (which is [K+]o dependent), γir is a term that includes the inward rectification properties of the channels, fATP is a term that expresses the ATP and ADP dependence of the current, Vm is membrane potential and EK the K+ equilibrium (Nernst) potential. The inward rectification term is expressed as follows:
where Kh,Mg and Kh,Mg are voltage and [K+]o dependent. Finally, the ATP-ADP dependence follows another Hill-type equation as follows:
where both the semi-inhibition constant KATP and the Hill exponent H depend on the intracellular ADP concentration.
2.2. Modeling drug-ion channel interaction
Also at the ionic channel level, the effect of antiarrhythmic drugs can be mathematically formulated. Models that simulate effects of drugs and reproduce the associated clinical electrocardiography signal are highly relevant to assess and predict antiarrhythmic actions, and help to improve the efficacy and safety of pharmacologic therapy. Additionally, it starts to be considered as an important technological tool for helping the pharmaceutical industry to develop new antiarrhythmic drugs spending less time and money. In this field, Dr. Starmer is one of the pioneers in modeling drug effects on sodium channels and an increasing number of research groups are developing 3D models of the heart to test drug actions. In the last years, our group has carried out different studies of antiarrhythmic and proarrhytmic effects of different drugs, such as pinacidil, dofetilide and lidocaine on cardiac arrhythmias.
To simulate the effects of drugs on ion channels, different kinds of models can be used (Brennan et al., 2009). Sometimes it is sufficient to change the conductance of ion channels by introducing the effect of the drug, but to be able to reproduce use- and frequency-dependent block it is necessary to use a state-dependent block model. In the state-dependent block model, drug interaction with ionic channels depends on the state of the channel. Drug-associated channels (blocked channels) do not conduct. The fraction of blocked non-conducting channels depends on the equilibrium between blocked and unblocked states. Two hypotheses have been proposed to reproduce the state-dependent block: the modulated receptor hypothesis and the guarded receptor hypothesis. The Modulated Receptor (MR) hypothesis (Hille, 1977; Hondeghem & Katzung, 1977) states that the affinity of a drug changes with the state of the channel. Additionally, the rate constants of state transitions are different for bound and unbound satates. This is the most general model to reproduce the effects of drugs; however a great number of parameters have to be estimated. The Guarded Receptor (GR) hypothesis (Starmer et al., 1984) assumes that the drug has a fixed and state-independent affinity for their receptor sites, but its access to the binding sites is controlled by the voltage-dependent channel gates. In this hypothesis, the rate constants for the state transitions do not change when the drug binds to the channel. The GR model is less complex than the MR model for the same drug-ion channel interaction as fewer parameters have to be estimated.
In the present section three different models of drugs: pinacidil, lidocaine, and dofetilide are described.
2.2.1 Modeling pinacidil effects
Pinacidil is a vasodilator drug and is known to increase IK(ATP) in cardiomyocytes, as well as in vascular smooth muscle and pancreatic Beta cells (Fan et al., 1990; Nakayama et al., 1990). The activation of ATP-sensitive potassium (KATP) channels in heart tissue, especially during myocardial ischemia, markedly changes action potential (AP) configuration by reducing AP duration (APD) (Nichols et al., 1991) and plays an important role in cardiac arrhythmogenesis (Janse & Kleber, 1981). To analyze the controversial effects of this drug through computer simulations, the first step is to build a model of the effects of pinacidil at the ion channel level. Pinacidil enhances the activity of IK(ATP) increasing the fraction of open KATP channels (fATP) (Fan et al., 1990). The equation of IK(ATP) formulated by Ferrero et al. (Ferrero, Jr. et al., 1996) was modified, as indicated in the following equations:
The equation of fATP was based on the pharmacokinetic scheme depicted in Fig. 1. This scheme describes the mechanism of action by which ATP molecules close KATP channels, while pinacidil molecules open it. ATP, ADP and pinacidil molecules bind to different sites in the protein receptor (Nichols et al., 1991; Weiss et al., 1992). According to the open states when molecules of pinacidil and/or ADP are bound to the protein channel, the fraction of activated KATP channels was expressed in equation (5). Experimental data (Fan et al., 1990; Nakayama et al., 1990; Weiss et al., 1992) were fitted to equation (5) to determine the values of equilibrium dissociation constants and the Hill coefficient (H) (see Trenor et al. 2005 for details). A factor (fPIN) addressing the voltage-dependent block of KATP channels by pinacidil was also included in IK(ATP) equation, as in Fan et al. (Fan et al., 1990).
2.2.2 Modeling lidocaine effects
Lidocaine is a class I antiarrhythmic drug that exerts its effect by blocking inward sodium current (INa) in a use-dependent manner, and is more effective for high stimulation frequencies (Clarkson et al., 1988). The block of lidocaine exists both as neutral and charged forms at a physiological pH. The charged form predominates at low pH (6.4) due to the protonation of the neutral form with hydrogen molecules. Experimental data have confirmed that lidocaine is more effective and recovery from block is slowed (Broughton et al., 1984; Moorman et al., 1986; Wendt et al., 1993) under acidosis.
The formulated model of lidocaine-INa interaction takes into account the modulatory effect of pH on the Na+ channel block and is based on the GR theory proposed by Starmer (Starmer & Courtney, 1986). Three main processes have been considered: the hydrophobic pathway, the hydrophilic pathway, and coupling between blocked channels by charged and neutral forms using a proton exchange process. Several experimental studies have suggested that the block of neutral drug takes place during the activated (O), inactivated (I) and resting (R) states (Bean et al., 1983; Clarkson et al., 1988), whereas the interaction between the charged form with the channel only takes place in the activated state (Liu et al., 2003; Moorman et al., 1986; Yeh & Tanguy, 1985). Fig. 2 represents the kinetic diagram of the complete model, and equations (6) and (7) mathematically describe the interaction of the neutral and charged form with Na+ channels, respectively.
bN and bC stands for the fraction of channels blocked by the neutral and charged form. The association and dissociation rate constants for the different states (activated, inactivated, and resting) are kO, kI, kR, lO, lI, and lR, respectively. The activation channel gate is m and the inactivaction channel gates are h and j. Additionally, we have taken into account the rate constants kp and lp to determine the proton exchange process. In equation (7), the association and dissociation rate constants are kC and lC, respectively.
As stated in the GR hypothesis, we introduced the factor (1-b), where b is the sum of the bN and bC in the INa equation.
2.2.3 Modeling dofetilide effects
Dofetilide is a specific and potent blocker of the rapid component of the delayed rectifier K+ current (IKr) with an IC50 in the nanomolar range (3.9-31 nmol/L) for ventricular myocytes (Carmeliet, 1992; Jurkiewicz & Sanguinetti, 1993; Weerapura et al., 2002b). Dofetilide is classified as a pure class III antiarrhythmic agent and provokes a prolongation of action potential duration (APD) without any effect on the resting membrane potential, action potential amplitude or maximum rate of depolarization (Jurkiewicz & Sanguinetti, 1993; Tande et al., 1990).
To model the effects of dofetilide on IKr, the formulation of IKr proposed by Zeng and Rudy (Zeng et al., 1995) for guinea pig ventricular cells was used. This formulation is based on experimental data obtained by Sanguinetti and Jurkiewicz (Sanguinetti & Jurkiewicz, 1990). In their model, IKr channels present three possible states: Closed (C), Open (O) and Inactivated (I). The formulation of IKr in the ventricular cell model (Zeng et al., 1995) is expressed in equation (8).
where Vm is the membrane potential, EKr is the reversal potential, GKrmax is the maximum conductance of IKr, Xr is the activation gate and R is the time-independent inactivation gate. In our dofetilide-IKr model, the effect of dofetilide is represented by introducing the factor (1b) in the IKr formulation (where b is the fraction of channels blocked by the drug). Thus, the new formulation of the rapid component of the delayed rectifier K+ current taking into account the effect of dofetilide is
The blocking activity (b) on IKr equation is modeled using the GR hypothesis (Starmer et al., 1984). Dofetilide experimental studies have shown drug-receptor interaction in open and inactivated states but not in closed states (Weerapura et al., 2002a; Yang et al., 1997). It has been also observed that the drug is trapped when the channel closes, thereby the transition between open and closed state is possible but dissociation of the drug does not occur at closed states (Carmeliet, 1992; Weerapura et al., 2002a). Taking into account these experimental evidences we suggest the model of IKr block by dofetilide shown in Fig. 3.
Using this model, the blocking factor b for a concentration of dofetilide D can be calculated as:
3. Modeling cardiac electrical cellular activity
The electrical activity of the cell is determined by the interaction of the multiple ionic channels and ion exchangers present in the cellular membrane. Their activity conditions the cardiac AP, which is intimately synchronized with the mechanical activity of the heart. Several computational models of cardiac cellular electrical activity have been developed, and characterize the activity of ventricular myocytes, atrial myocytes, sinus node and Purkinje fibers for different animal species. The degree of complexity of these models increases going from heuristic models describing the physical process of cell excitation and propagation, such as Fitz-Hugh-Nagumo model (Fitzhugh, 1961), to phenomenological models (Bueno (Bueno-Orovio et al., 2008) model), and to extremely detailed ionic models. The most detailed ionic models for AP include transmembrane currents, ion transfer through transporters and dynamic changes of ion concentrations during the action potential. The ventricular AP models most extensively used in simulations because of their detailed and realistic formulation are Luo and Rudy model (Faber & Rudy, 2000; Luo & Rudy, 1994; Zeng et al., 1995) for guinea pig, Shannon (Shannon et al., 2005) and Mahajan (Mahajan et al., 2008) models for rabbit, and ten Tusscher (ten Tusscher & Panfilov, 2006) and Grandi (Grandi et al., 2010) models for human. As regards atrial myocytes, Nygren (Nygren et al., 1998) model and Courtemanche (Courtemanche et al., 1998) model have also had an important repercussion. Luo Rudy AP model was one of the first detailed models formulated and has been a fundamental base on which other models have been built. This model has been improved and extensively used to simulate the ventricular electrical activity. In these detailed AP models, the formulation of the ionic currents follow Hodgkin-Huxley, Goldmann-Hodgkin-Katz (GHK) or Markov formalisms, and the membrane potencial Vm is related to the ionic currents through equation (11).
where the total ionic current through the cell membrane (Im) is expressed as the sum of all ionic currents (Is through ionic channels and Ip through pumps and exchangers) plus the current through the membrane equivalent capacitor.
Fig. 4 shows a scheme of the main ionic currents of a ventricular myocyte (panel A) and the time course of membrane potential and some ionic currents (panel B) obtained by solving the non-linear differential equations system defined by equation (11) and the different ionic current equations in Luo and Rudy model (Faber & Rudy, 2000; Zeng et al., 1995).
3.1. Modeling myocardial ischemia at cellular level
Acutely ischemic action potentials can be simulated using these AP models appropriately modified to represent ischemia-activated currents (such as IK(ATP)) and ischemic-modified currents (e.g. the inward Na+ and Ca2+ currents), as explained previously.One of the main components of acute ischemia is hypoxia (i.e. lack of oxygen in the tissue), and the main change exerted by hypoxia to myocardial action potentials is a significant shortening in APD (Morena et al., 1980). The model proposed by Ferrero et al. for the ATP-sensitive K+ current (Ferrero, Jr. et al., 1996) yielded the first theoretical proof of the determinant role of the IK(ATP) current in hypoxic action potential shortening (Ferrero, Jr. et al., 1996), which is an arrhythmic factor due to the APD dispersion it generates across the myocardium (Janse & Wit, 1989). Indeed, as shown in Fig. 5, although very few KATP channels do activate during acute myocardial ischemia (<1% in the first 15minutes), this degree of activation is enough to account for the >50% shortening in APD which is observed in experimental conditions (Morena et al., 1980).
When the simulated myocyte is subject to hyperkalemia and acidosis, as well as hypoxia, the main features of the ischemic AP waveforms are correctly reproduced (Ferrero, Jr. et al., 2003; Shaw & Rudy, 1997a). A shown in Fig. 6, ischemia produces APD shortening, resting potential elevation, upstroke phase division into two phases (the first one mediated by inward Na+ current, the second one by inward Ca2+ current through the L-type channels). Moreover, acute ischemia changes important intrinsic AP properties which are pivotal in arrhythmogenesis (e.g., postrepolarization refractoriness). Overall, this ischemic AP simulations are the first level of the multi-scale ischemia simulations that will be described further in this chapter.
3.2. Effects of drugs at cellular level
As mentioned before, the electrical activity of the cell is altered under the effects of drugs. Once the effects of the drug have been formulated at the ionic channel level, the alterations of the electrical activity of the cell can be analyzed including the modified current in the AP model. In this way, the effects of pinacidil and lidocaine on ischemic cells are described. Also the effects of dofetilide on different types of ventricular cells will be analyzed.
3.2.1. Effects of pinacidil on AP
As mentioned in section 2, pinacidil shortens the APD by the activation of KATP channels. To validate pinacidil model, the formulation of the effect of pinacidil (equations (4) and (5)) on IK(ATP) was introduced into the Luo and Rudy model (Luo & Rudy, 1994; Zeng et al., 1995). Simulations were carried out combining different concentrations of pinacidil with various [ATP]i and [ADP]i values corresponding to several stages of early acute ischemia (Cole et al., 1991; Nakayama et al., 1990; Weiss et al., 1992). The resulting APs are shown in Fig. 7. Consistent with experimental findings (Cole et al., 1991; Nakayama et al., 1990), our results show that increasing concentrations of pinacidil and/or decreasing intracellular ATP levels reduce APD. It is noteworthy that changes concerning intracellular ATP and ADP concentrations, occurring during acute ischemia, significantly enhance the effects of pinacidil. These changes in APD (and hence in cell refractoriness) could be key in determining the reentry probability in acute ischemia after the administration of the drug.
3.2.2. Effects of lidocaine on AP
The block of INa by lidocaine leads to a decrease in the maximal AP upstroke (dV/dtmax). To simulate the effects of lidocaine on AP features, we incorporated the model of lidocaine into the Luo and Rudy AP model (Luo & Rudy, 1994). Specifically, we measured the reduction of the maximum AP upstroke for different stimulation frequencies and pH values when increasing concentrations of the drug. We examined the response of dV/dtmax to changes in pH from 7.4 to 6.4 for different concentrations of the drug (50 and 100 of µmol/L lidocaine) and at different stimulation frequencies. Panel A in Fig. 8, shows a decrease of 66% and 83% in dV/dtmax from the control value for 50 and 100 µmol/L lidocaine, respectively, at a basic cycle length (BCL) of 500 ms and pH 6.4. For a pH of 7.4 (see Fig. 8 panel B), dV/dtmax was reduced in 45% and 61% for the same concentrations and BCL values. These data support the fact that lidocaine has a greater effect under lower pH. Furthermore, the use-dependent effect is more pronounced at lower pH, as the differences in the decrease of dV/dtmax for different stimulation frequencies are higher for low pH.
3.2.3. Effects of dofetilide on AP
The steady-state effects of dofetilide concentration on AP characteristics for the different kinds of guinea pig ventricular cells at different BCL values are shown in Fig. 9. In all simulations, once dofetilide was applied the cellular models were stimulated until steady- state was reached. Only magnitudes associated to the last stimulation pulse are shown. Fig. 9A shows that the prolongation of APD90, for all tested dofetilide concentrations, significantly depends on the type of cell. APD90 of M-cells was markedly prolonged compared to endocardial and epicardial cells. For a concentration of 10 nmol/L, APD90 of endocardial cells was prolonged in 8.8%, 11.5% and 12.4% of control value (APD90 before the application of drug) for BCL of 300, 1000 and 2000 ms, respectively. For the same conditions, APD90 of epicardial cells was increased in 5.8%, 6.9% and 7.1% of control, respectively. However, under these conditions, APD90 of M-cells was substantially increased in 14.8%, 36.1% and 292.7% of control, respectively. The three types of cells present reverse use-dependency. Due to this effect, the prolongation of APD90 increased whith the BCL. Fig. 9B shows the relative increase of APD90 (ΔAPD90 in%) induced by different concentrations of dofetilide at different BCLs.
Fig. 9B shows that dofetilide presents reverse use-dependence in the range of tested concentrations (10 nmol/L to 1 µmol/L) with a much more pronounced effect in M-cells than endocardial or epicardial cells. Even more, when the concentration of dofetilide increased, the reverse use-dependency was enhanced.
4. Modeling cardiac tissue electrical activity
Abnormalities in impulse propagation are involved in the genesis and maintenance of cardiac arrhythmias and many experimental and theorical studies have analyzed action potential propagation both in normal and pathological situations. The electrical coupling between cells and the propagation process is described in the bidomain model (Tung, L., 1978). This model consists of two continuous domains, the intracellular and the extracellular domain, where electrical currents are governed by Ohm´s law. The following equations describe the bidomain model:
where Di and De are the intracellular and extracellular conductivity tensors, respectively, Ve and Vm stand for the extracellular and the membrane potentials, respectively, Cm is membrane capacitance, Iion represents the total ionic current density computed using AP models, and Ist is the stimulus current density.
Considering that anisotropy similarly affects the intracellular and the extracellular conductivity tensors (Di=De, being a scalar variable), the bidomain model can be simplified into the monodomain model and the previous equations can be reduced to the following one:
Subsequently, the monodomain model consists of an elliptic partial differential equation and a parabolic partial differential equation coupled to a system nonlinear ordinary differential equations describing the ionic current through the cellular membrane. As this model is mathematically simpler and less computational demanding than the bidomain model, it is widely used for cardiac electrophysiology simulations. Heterogeneous structures of cardiac tissue can also be considered in tissue models allowing a wide use of cardiac models in the study of ionic mechanims of arrhythmias.
4.1. Tissue model of myocardial ischemia
In the case of pathological tissues, not only the altered AP models have to be considered but also the spatial inhomogeneities in the electrophysiological changes distributed within the tissue. If we consider the ischemic pathology, after coronary occlusion, the lack of oxygen and blood flow provokes important and heterogeneous electrophysiological changes in the affected ventricular cells defining an ischemic zone and a border zone, setting the stage for reentrant arrhythmias and VF. A large body of clinical and experimental studies try to unravel the mechanisms by which acute myocardial ischemia provokes life-threatening arrhythmias. However, the rapid changes arising during acute myocardial ischemia and the limitations in experimental techniques hamper the complete understanding of this pathology. In this way, computer simulations overcome these limitations providing a helpful tool to understand the electrophysiological changes arising during myocardial ischemia and the underlying mechanisms of reentrant arrhythmias observed during its acute phase. At the tissue level, our group has developed a regional model of ischemia, and the consideration of the temporal evolution of ischemia has allowed undertaking computational studies to assess the changes in the vulnerability of the ischemic tissue to reentry. Not only ionic currents have been analyzed to understand the process of reentry initiation but also the traditional indicator for reentry, i.e. dispersion of repolarization, and the safety factor (SF) for conduction (a parameter that quantifies the source-sink relationship of propagation) as an indicator for conduction block, have been evaluated.
A regional phase Ia ischemic 2D anisotropic monodomain tissue (Ferrero, Jr. et al., 2003; Romero et al., 2009; Trenor et al., 2007) was modeled with a high degree of electrophysiological detail. Changes in electrophysiological parameters were modeled according to experimental results (Coronel, 1994; Ferrero, Jr. et al., 1996; Irisawa & Sato, 1986; Weiss et al., 1992; Yatani et al., 1984). These changes are described in Fig. 10A, which shows the time-course of the three main components of acute ischemia: hypoxia, hyperkalemia and acidosis. Fig. 10B represents the distribution of a simulated cardiac 2D tissue comprising a normal zone (NZ), a circular central ischemic zone (CZ), and a ring-shaped ischemic border zone (BZ) surrounding the CZ. In the BZ all parameters affected by ischemia follow linear spatial gradients, as experimentally recorded (Coronel et al., 1988). Membrane kinetics was simulated by using a modified version of the Luo-Rudy dynamic AP model (LRd00) (Faber & Rudy, 2000) (see Romero et al. 2009 for details), which includes the mathematical formulation of IK(ATP) proposed by Ferrero et al. 1996 (Ferrero, Jr. et al., 1996).
A basic stimulus (S1) followed by a prematures stimulus (S2) at different coupling intervals (CIs) was applied on the edge of the tissue to quantify the vulnerable window (VW), an indicator of the vulnerability to reentry. The VW was defined as the interval of CIs that led to reentry. Fig. 10C shows gray-coded voltage snapshots of a figure-of-eight reentry developed after premature stimulation at minute 8.25 of ischemia. Our results demonstrated that during myocardial ischemia phase Ia, the VW of a bidimensional tissue had an asymmetric unimodal distribution, peaking at minute eight after coronary occlusion (see Fig. 10D), which is in close agreement with experimental observations (Cascio, 2001). We also analyzed the SF for conduction in the regional ischemic tissue using an improved version (Romero et al., 2009) of the formulation proposed by Shaw and Rudy (Shaw & Rudy, 1997c). Our results also indicate that, as ischemia progresses, the SF decreased and the heterogeneity of refractoriness grows. Therefore, a lack of correlation between dispersion of refractoriness and vulnerability to reentry was observed, although our results corroborate the fact that dispersion of refractoriness is essential to reentry generation. Moreover, we found that, in approximately half of the simulated reentries, the area where unidirectional block (UDB) took place was completely recovered from refractoriness when the AP propagation failed. However, the line where the SF dropped below unity matched very much the area where propagation failed in every simulation. Thus, it seemed that the reduction of the source-sink ratio is the ultimate cause of the UDB leading to reentry.Reentrant patterns of activation during regional acute ischemia have been mechanistically analyzed also under the effects of pharmacologic agents, such as pinacidil or lidocaine.
4.2. Drugs effects on reentry during myocardial ischemia
To analyze the arrhythmogenic effects of pinacidil, different sets of simulations were carried out using different concentrations of the drug. In each set the regionally ischemic 2D tissue described above was stimulated using the aforementioned S1-S2 protocol and the VW was evaluated. Table 1 shows that increasing the concentration of pinacidil had an interesting non-monotonic effect on the duration of the VW. Indeed, low concentrations of the drug increased the VW, until a maximum of 54 ms was reached for [P]=3 m/L. If the concentration was further increased, the VW decreased, until the VW vanished for [P]=10 m/L.
The simulations allowed to unravel the mechanisms underlying pinacidil-induced changes in vulnerability, through the analysis of the patterns of excitation in the tissue and the time evolution of ionic currents and membrane potential. Our results showed the existence of two opposite effects of pinacidil in terms of AP propagation, which modulate the biphasic effects on the VW. On the one hand, pinacidil enhances IK(ATP) (Ferrero, Jr. et al., 1996), which counteracts inward currents responsible for depolarization, mainly INa (which is depressed as a consequence of ischemia) and ICa(L), so that AP reaches lower values of Vm during the depolarization phase. In this way, the potential gradient in the direction of propagation is also reduced, and less electrotonic current flows to downstream cells (Shaw & Rudy, 1997a), reducing axial current. In these conditions, the SF is reduced (Shaw & Rudy, 1997c) and propagation failure may develop. On the other hand, the IK(ATP) enhancement provoked by pinacidil reduces APD, which results in an earlier recovery of excitability (Shaw & Rudy, 1997b), favoring AP propagation. Both effects of pinacidil are related as follows: the earlier the cell recovers its excitability, the less axial current is needed to elicit an AP.
We also simulated the effects of lidocaine on a regionally ischemic bidimensional tissue prone to reentrant circuits. Simulations were also carried out using an S1-S2 protocol applied on one edge of the tissue and the VW for reentry was calculated for various degrees of severity of ischemia and various concentrations of the drug (20, 50, and 100 µmol/L). The VWs for reentry are shown in Fig. 11, and follow an unimodal behavior with ischemia time course, peaking (35 ms) for 6 minutes after the onset of ischemia. When the concentration of lidocaine was increased, the VW became wider, indicating a higher vulnerability to reentry in the presence of the drug. This effect can be due in part to the heterogeneous action of the drug in the diverse zones of the tissue, contributing to the dispersion of CV and setting the stage for reentry.
4.3. Effects of dofetilide on 1D heterogeneous tissue
In the case of dofetilide, a class III antiarrhythmic drug recently approved by FDA (Food and Drug Administration) for the treatment of persistent atrial fibrillation and flutter, different studies have questioned the antiarrhythmic action of dofetilide in preventing and terminating ventricular tachycardias. Indeed, dofetilide promotes the prolongation of QT interval, which has been related to the trigger of a polymorphic ventricular tachycardia called torsade de pointes. The study of the proarrhythmic effects of dofetilide has also been undertaken with computer simulations at the tissue level, where pseudo ECGs have been evaluated, giving understanding about the mechanisms by which QT interval is altered by this drug.
To study the effect of dofetilide on transmural dispersion of repolarization (TDR), a linear strand including the three kinds of cells was used. Fig. 12A illustrates the steady-state APD90 distribution along a well-coupled heterogeneous fiber and the steady-state APs time-course of the three types of cells belonging to the different regions of the tissue (right side in Fig. 12) under two different conditions: before (control) and after the application of dofetilide 100 nmol/L. The pseudo ECGs are also represented before and after the application of the drug.Our results show that in a well-coupled strand under control conditions, at a BCL of 1000 ms, M-cells present a maximum APD90 of 191 ms, which is 15 ms and 35 ms longer than the shortest APD of endocardial (176 ms) and epicardial (156 ms) cells, respectively. The application of dofetilide 100 nmol/L increased APD90 dispersion (difference between the maximum and the minimum APD90 in the fiber) from 35 ms to 58 ms, being the maximum APD90 recorded in the M-zone and the minimum in the epicardial zone. The electrograms show an increment in the QT interval induced by the application of dofetilide from 220.8 ms (control) to 274.6 ms for 100 nmol/L dofetilide.
The effect of decreasing the intercellular coupling (intercellular coupling resistance Ri multiplied by three) at a BCL of 1000 ms is shown in Fig. 12B. Poor coupling increased the APD90 dispersion from 35 ms (normal coupling) to 58 ms. Even more, under poor coupling conditions dofetilide induced a steeper distribution of APD along the fiber and the APD90 dispersion increased to 102 ms, for dofetilide concentration of 100 nmol/L. The QT intervals were higher than the ones observed in a well coupled fiber. The application of 100 nmol/L dofetilide induced a sharper increase of QT from a value of 251.8 ms in control, to 318 ms.
5. Modeling whole heart electrical activity
Whole heart models are needed to study arrhythmias that critically depend on the spatial organization of the heart. Some spatial physiological features may be localization dependent, such as regional ischemia or localized infarct scars. Furthermore, certain arrhythmias are uniquely encountered in the whole heart, such as large reentrant circuits in acute ischemia or initiation of atrial fibrillation by pulmonary venous foci. Other important considerations of the whole heart models are the cardiac fibers arrangement and transmural inhomogeneity in repolarization. Cardiac fibers are arranged as counter-wound helices encircling the ventricular cavities, and the local orientation of these fibers depends on transmural localization (Helm et al., 2005). On the other hand, the heterogeneous expression of IKs, Ito, and INCX, among others, leads to transmural inhomogeneity in repolarization (Antzelevitch & Fish, 2001; Antzelevitch et al., 1991). The vertiginous progress in medical images allows building realistic models of the heart including its highly complex anatomical structure of ventricles and/or atria (see Clayton et al. 2011 for review).
The mathematical problem associated to the resolution of the differential equations modeling the electrical propagation in a 3D heart model does not have an analytical solution. Indeed, electric conduction is described by a set of partial differential equations (PDEs) and ionic currents through the cell mambrane are described with a nonlinear stiff system of ordinary differential equations (ODEs). The numerical solution of the equations is very computationally demanding (Heidenreich et al., 2010b) and require the discretization in space and time of PDEs, as well as the integration of nonlinear systems of ODEs. The operator splitting technique for the time discretization (Qu & Garfinkel, 1999) is used to solve the equations system. The numerical methods used, such as finite differences method or finite elements method, are computationally demanding and require the use of high performance computing techniques.
Whole heart models also permit to relate the arrhythmic behavior with its manifestation in electrograms on the cardiac surface and in electrograms on the body surface, through a torso model (Weiss et al., 2009).
5.1. 3D simulations of myocardial ischemia
The use of these ventricular 3D models has shed light into the mechanisms by which reentrant arrhythmias are initiated during regional myocardial ischemia. As an example, a 3D model of the human ventricles has been recently used to study the appearance of figure-of-eight reentry under regional acute ischemic conditions (Heidenreich et al., 2010a; Heidenreich et al., 2010b).
Fig. 13 shows different views of the simulated human ventricles. The detailed realistic local fiber orientation obtained from DTI images (diffusion tensor imaging) was introduced in the model and is key in determining the electrical propagation. An ischemic region was introduced in the left ventricle, following the same normal-border-central zone scheme explained in a previous section of this chapter. The central and right panels in Fig. 13 depict the values and gradients of the main ischemic parameters (namely [K+]o, [ATP]i, INa maximum conductance and ICa(L) maximum conductance).
As shown in Fig. 14, the model predicted the generation of figure-of-eight reentries, which cross the central ischemic zone formed in the epicardial surface due to the longer refractory period of the midmyocardial layers. Also, focal activity experimentally observed in the epicardium could be caused by re-entrant wavefronts propagating in the mid-myocardium that reemerge in the heart surface.
Multiscale modeling of the electrical cardiac activity represents nowadays a very valuable tool in cardiac disease research which complements experimental and clinical studies. Although significant progresses are made in this field, further improvements are required, such as the modeling of electromechanical coupling of the heart activity, which currently represents an important challenge for modelers. The vertiginous technical improvements in the obtention of medical images and genetic, molecular and ionic measurements, represent a great help to the modeling field, which is in the process of creating a virtual human.
This work was partially supported by the European Commission preDiCT grant (DG-INFSO-224381), by the Plan Nacional de Investigación Científica, Desarrollo e Innovación Tecnológica del Ministerio de Ciencia e Innovación of Spain (TEC2008-02090; TIN2004-03602; TEC 2005-04199/TCM; TIC 2001-2686), by the Programa de Apoyo a la Investigación y Desarrollo (PAID-06-09-2843) de la Universidad Politécnica de Valencia, by the Dirección General de Política Científica de la Generalitat Valenciana (GV/2010/078).