Mathematical Model of an Adult Human Atrial Cell
The Role of K+ Currents in Repolarization
Abstract—We have developed a mathematical model of the human atrial myocyte based on averaged voltage-clamp data recorded from isolated single myocytes. Our model consists of a Hodgkin-Huxley–type equivalent circuit for the sarcolemma, coupled with a fluid compartment model, which accounts for changes in ionic concentrations in the cytoplasm as well as in the sarcoplasmic reticulum. This formulation can reconstruct action potential data that are representative of recordings from a majority of human atrial cells in our laboratory and therefore provides a biophysically based account of the underlying ionic currents. This work is based in part on a previous model of the rabbit atrial myocyte published by our group and was motivated by differences in some of the repolarizing currents between human and rabbit atrium. We have therefore given particular attention to the sustained outward K+ current (Isus), which putatively has a prominent role in determining the duration of the human atrial action potential. Our results demonstrate that the action potential shape during the peak and plateau phases is determined primarily by transient outward K+ current, Isus, and L-type Ca2+ current (ICa,L) and that the role of Isus in the human atrial action potential can be modulated by the baseline sizes of ICa,L, Isus, and the rapid delayed rectifier K+ current. As a result, our simulations suggest that the functional role of Isus can depend on the physiological/disease state of the cell.
Within the past 5 years, quite extensive voltage-clamp and action potential data from human atrial myocytes have been published from a number of laboratories. In most cases the human atrial tissue was obtained during open heart surgery, in which small pieces of the right atrial appendage are excised as part of the cannulation procedure for cardiopulmonary bypass. However, cardiac tissue can also be obtained from the free wall of the atrium during valve replacement procedures, from explanted failing hearts during heart transplantations, or from donor hearts that cannot be used for transplantation. Once a specimen of cardiac tissue is obtained, enzymatic dispersion techniques are used to isolate single cardiac cells, from which action potential and voltage-clamp data can be recorded. We have made recordings of the outward K+ currents, which are responsible for repolarization in human atrial myocytes.1 2 On the basis of these data, as well as other published results, we have developed the first comprehensive mathematical model of the electrophysiological responses of a representative human atrial cell.
There are a number of published mathematical models that simulate the electrophysiological responses in several different species and cardiac cell types. Examples include the Purkinje fiber model of DiFrancesco and Noble,3 the Hilgemann and Noble atrial model,4 the Earm and Noble model of the single atrial cell,5 the bullfrog atrial and sinus venosus models of Rasmusson and colleagues,6 7 the ventricular cell models of Luo and Rudy,8 9 10 11 the rabbit sinoatrial node cell model of Demir et al,12 and the rabbit atrial cell model of Lindblad et al.13 Recently, the emphasis has shifted from general models, based on voltage-clamp data from several species,3 4 to more detailed models based on data from single isolated cells from a particular species. This is a direct reflection of the progress made in experimental work and the resulting availability of more comprehensive data. Our goal was to develop a model that is sufficiently accurate to have predictive capabilities for selected aspects of the electrophysiological responses in human atrium. Emphasis has been placed on the functional roles of the K+ currents during repolarization.
Materials and Methods
Our model is of the general type first introduced by DiFrancesco and Noble3 and is based, in part, on the rabbit atrial model of Lindblad, Murphey, Clark, and Giles13 (hereafter referred to as the LMCG model). As shown in Fig 1⇓, the model consists of a Hodgkin-Huxley–type electrical equivalent circuit for the sarcolemma coupled with a fluid compartment model. The dimensions of the human atrial myocyte are assumed identical to those of the rabbit atrial cell in the LMCG model. These dimensions (cylindrical geometry of 130-μm length and 11-μm diameter) are very close to the dimensions of human atrial myocytes (eg, 120-μm length and 10- to 15-μm diameter14 ). In addition, we have used a total cell capacitance of 50 pF,which agrees very well with our experimental observations for human atrial myocytes (51.9±3.5 pF, n=52).
Fig 1A⇑ shows the electrical equivalent circuit for the sarcolemma of the human atrial cell. It includes each of the ionic currents that are known to contribute to the action potential in human atrial myocytes (INa, ICa,L, It, Isus, IK,r, IK,s, and IK1), the Ca2+ and Na+-K+ pump and Na+-Ca2+ exchanger currents responsible for maintaining intracellular ion concentrations (ICaP, INaK, and INaCa), and the Na+ and Ca2+ background (leakage) currents (IB,Na and IB,Ca). Mathematical expressions describing the time and voltage dependence of the ionic currents have been developed on the basis of published voltage-clamp data recorded predominantly from human atrial myocytes. (See “Glossary” after Appendix 2 for terms used in text, figures, and tables.)
Na+ Current: INa
Voltage-clamp data for INa in human atrial and ventricular myocytes, recorded at 17°C have been published by Sakakibara and colleagues.15 16 These data suggest that the activation threshold is very close to the resting potential (≈−75 mV), which seems unrealistic, given that atrial and ventricular cells exhibit stable resting potentials and thresholds for activation near −55 mV. Moreover, the steady-state inactivation curves measured by Sakakibara and colleagues15 16 are such that INa would be completely inactivated (ie, no current available) at the resting potential. As pointed out by these authors, this is probably a result of time- and/or temperature-dependent shifts in the steady-state inactivation characteristics. Other results, such as the data from rabbit atrium published by Wendt et al17 (on which the INa description in the LMCG model is based) yield more positive (depolarized) steady-state activation and inactivation curves.
In developing a model of INa under physiological conditions, we have found it necessary to use indirect information about INa provided by action potential data in addition to voltage-clamp data. Thus, we have adjusted the steady-state activation curve (m̄3) for INa so that the threshold at which an action potential is elicited agrees with experimental observations. Furthermore, the peak magnitude of INa was adjusted to match the maximum upstroke velocity of the action potential. (For a discussion of the relation between INa and action potential upstroke velocity, refer to Cohen and Strichartz.18 ) Fig 2A⇓ shows the steady-state activation (m̄3) and inactivation (h̄) curves used to model INa. Compared with the data from human atrial myocytes obtained by Sakakibara et al15 (Fig 2A⇓), there are significant positive shifts in both the m̄3 (+22.8 mV) and the h̄ (+32.2 mV) curves. Fig 2B⇓ shows a simulated peak current-voltage relationship (steps from a holding potential of −80 mV) for INa. The mathematical expressions for the kinetics of activation (τm) and inactivation (τh1 and τh2) are very similar to those of the LMCG model, as shown in Figs 2C⇓ and 2D⇓. The processes of inactivation and recovery from inactivation are both described by the sum of a fast and a slow exponential. At plateau potentials, the fast component of inactivation has a time constant of 0.3 ms. The slow component of inactivation accounts for 10% of the total current and has a time constant of 3.0 ms at plateau potentials; ie, it is 10 times slower than the fast component, in agreement with the data of Sakakibara et al.15
L-Type Ca2+ Current: ICa,L
Several quite comprehensive studies of ICa,L in human atrial myocytes have been published.19 20 21 22 23 24 25 Overall, the results are quite consistent: ICa,L has an activation threshold near −40 mV, its peak at ≈0 mV, and an apparent reversal potential of +50 to +60 mV. The current density for ICa,L varies considerably among these studies, however, and is also known to be reduced in diseased or dilated cells.22 26 Table 1⇓ shows published ICa,L densities compared with those used in the present model. The inactivation process for ICa,L is usually described as the sum of a fast and a slow component.19 20 21 22 23 In order to accurately measure the steady-state voltage dependence of inactivation, the duration of the “conditioning” prepulse in the voltage-clamp protocol (eg, see Li and Nattel25 ) must be at least four to five times the time constant of the slower component of inactivation. If the prepulse duration is shorter, the slower component will not have reached its steady state when the availability of current is measured (second voltage-clamp pulse), and inactivation will appear incomplete. Ouadid et al24 measured inactivation at room temperature (where the slower time constant exceeds 100 ms23) and found that inactivation was incomplete (U-shaped inactivation curve) when the prepulse duration was 150 ms but that it became more complete if the prepulse duration was increased to 400 ms. The data obtained by Li and Nattel at 36°C yield almost complete inactivation characteristics for a prepulse duration of 150 ms, consistent with faster kinetics at this temperature. However, these authors also note that prolonging the prepulse results in more complete inactivation curves. On the basis of these results, we have modeled inactivation of ICa,L as a process involving two components with different time constants but identical, fully inactivating, steady-state voltage dependence. This formulation differs from the LMCG model13 as well as the Luo-Rudy model,9 which both use a single (voltage-dependent) component of inactivation with incomplete steady-state voltage dependence.
Ouadid et al24 also report that the inactivation process is slowed considerably when Ba2+ is substituted for Ca2+ as charge carrier, indicating that inactivation is also Ca2+ dependent. This phenomenon is included in the Luo-Rudy model9 but not in the LMCG model.13 Our formulation differs from that of Luo and Rudy, however, in order to incorporate recent experimental results. First, there is now evidence that there exists a small restricted subsarcolemmal domain between the L-type Ca2+ channels and the peripheral junctional sarcoplasmic reticulum (SR) where Ca2+ concentration ([Ca2+]d) may transiently reach much higher levels than in the cytosol as a whole.27 This subsarcolemmal domain is not included in the Luo-Rudy model, which models Ca2+-dependent inactivation as a function of total cytosolic Ca2+ concentration. We have included such a subsarcolemmal domain and have modeled Ca2+-dependent inactivation as a function of [Ca2+]d. Second, there is recent evidence to suggest that [Ca2+]d modulates ICa,L inactivation by promoting a rapid mode of inactivation.28 29 Therefore, in our model, [Ca2+]d determines the fraction of L-type channels that are in the rapidly inactivating mode, ie, the ratio of the fast to slow components of inactivation discussed above. Thus, in this model, an experiment with Ba2+ as charge carrier shifted this equilibrium so that all L-type channels were in the slow mode. With Ca2+ as charge carrier, however, the equilibrium is shifted toward the faster mode of inactivation by an amount determined by [Ca2+]d. Hence, the inactivation of the total ICa,L follows a biexponential time course, where the relative contributions of the fast and slow exponentials are determined by [Ca2+]d. We have chosen to model this [Ca2+]d dependence as an instantaneous function of [Ca2+]d, assuming that the shift between the two modes is rapid compared with the diffusion of Ca2+ out of the restricted domain.
Fig 3A⇓ shows the steady-state activation (dL) and inactivation (fL) curves used to model ICa,L. Assuming that activation can be measured more accurately at room temperature than at physiological temperature, we have used the activation curve (dL) measured by Mewes and Ravens.23 However, in order to fit the voltage-clamp data (peak currents) of Li and Nattel25 (Fig 3B⇓), we found it necessary to shift this activation curve by +3 mV. The expression for the inactivation curve (fL) is identical to that reported by Li and Nattel for human atrial myocytes at physiological temperature. Furthermore, the reversal potential for ICa,L is set to a constant value of +60.0 mV, as measured by Li and Nattel, rather than to the Nernst potential for Ca2+ ions. As shown in Fig 3B⇓, the simulated peak current-voltage relationship (steps from a holding potential of −80 mV) agrees well with voltage-clamp data. Figs 3C⇓ and 3D⇓ show the time constants of activation (τdL) and inactivation (τfL1 and τfL2) plotted against membrane voltage. We have used the inactivation and recovery time constant data obtained by Li and Nattel to formulate the expressions for τfL1 and τfL2 (see Fig 3D⇓). These data were obtained in human atrial myocytes at physiological temperature. The expression for τdL (Fig 3C⇓) is similar to that of the LMCG model.
Transient and Sustained Outward K+ Currents: It and Isus
Voltage-clamp experiments designed to study the outward currents responsible for repolarization in human atrial myocytes have identified a transient outward K+ current (denoted It), which activates rapidly on depolarization.1 30 In addition to this (Ca2+-independent) transient K+ current, a Ca2+-dependent transient outward current, which is activated by relatively large increases in [Ca2+]i, is sometimes observed.31 We have chosen not to include this Ca2+-dependent current, since it has never been observed in our experimental work. After It has decayed (inactivated), a more slowly inactivating, “sustained,” outward K+ current (denoted Isus) is observed.2 22 32 33 The available data suggest that Isus is a separate current from It and that it is also carried mainly by K+ ions.2 34 35
The literature regarding the voltage dependence of It is somewhat inconsistent. In particular, the results regarding steady-state activation vary considerably. Shibata et al1 report a half-activation voltage (V1/2) of +1.0 mV for It in human atrium, whereas Näbauer et al,36 Wettwer et al33 (both in studies of human ventricle), and Le Grand et al22 (human atrium) report values of +16.7, +20.6, and +33.3 mV, respectively. These results seem unrealistic, since It is known to have a strong influence on early repolarization,1 2 and these values of V1/2 would result in only a small amount of It current being activated during a normal action potential (peak at +20 to +30 mV). Among the possible explanations for this variability in V1/2 are the sensitivity of this parameter to the concentration of divalent cations (Cd2+ and Co2+) used to block ICa,L33 37 and the possibility of a difference between human atrial and ventricular It. Therefore, we have based our model on the data of Shibata et al,1 which were recorded from human atrium in the presence of a low (100 μmol/L) [Cd2+]. The data regarding the voltage dependence of steady-state inactivation are also variable between different studies. Fig 4A⇓ shows the steady-state activation (r̄) and inactivation (s̄) curves used to model It. The r̄ curve is based on a fit to data from Shibata et al (○ in Fig 4A⇓), and the s̄ curve is that reported by Firek and Giles.2 Time constants of activation (τr) and inactivation (τs) are plotted against voltage in Fig 4C⇓ and 4D⇓, respectively. Data regarding the time constant of inactivation2 indicate that τs is ≈13 ms (at 33°C) at membrane voltages positive to 0 mV. At hyperpolarized potentials, the recovery of It from inactivation as measured in our laboratory appears to follow an exponential time course with a strongly voltage-dependent time constant, increasing from 15 ms at −100 mV to 387 ms at −60 mV. Our formulation for τs is based on a fit to experimentally obtained recovery time constant values at negative potentials and inactivation time constant values at positive potentials. The recovery of It from inactivation in human atrial cells has been shown to be considerably more rapid than in rabbit atrial cells.2 38 As a result, It magnitude and action potential waveshape are much less rate dependent in human atrial cells.
Our model of Isus is based on the data of Wang et al.35 Fig 4B⇑ and 4D⇑ show the steady-state activation curve (rsus) and time constant of activation (τrsus) for this current. On the basis of our recent data, we have also included a slow (τ≈300.0 ms) partial (40%) inactivation. The time constant for this inactivation process was obtained by fitting a biexponential function to the decaying outward current waveforms as described by Koidl et al.39 In addition to providing an estimate of τssus, this method also provides a more accurate separation of It and Isus than when Isus is estimated as the current at the end of the pulse.2 35 40 By combining our models of It and Isus, we are able to produce current waveforms closely resembling those recorded from human atrial myocytes in response to voltage-clamp pulses. Fig 5⇓ shows how It (panel A) and Isus (panel B) combine (panel C) to produce waveforms similar to experimental results (panel D). In our experience, the size of these currents varies considerably between individual cells, and as reported by Amos et al,40 there is also considerable variability in the ratio of It to Isus. The sizes of It and Isus in the model were chosen to provide good fits to action potential data and are well within this experimental variability. Table 2⇓ compares It and Isus densities reported in the literature with those used in this model.
Delayed Rectifier K+ Currents: IK,r and IK,s
Recent studies of the delayed rectifier K+ currents in human and rabbit atrial myocytes show that in both species the delayed rectifier current is generated by two distinct K+ conductances. These “rapid” (IK,r) and “slow” (IK,s) conductances have significantly different properties and can be separated experimentally on the basis of, for example, the sensitivity of IK,r, but not IK,s, to the antiarrhythmic drug E-4031.41 42 43 Since IK,s is believed to contribute only a small fraction of the total delayed rectifier current during a normal atrial action potential,41 one could produce an acceptable fit to nominal action potential data using a model incorporating only IK,r. However, because of the very slow activation characteristics of IK,s, this current would be expected to be more significant at high heart rates, where it could build up progressively from cycle to cycle as a result of residual activation, ie, failure to decay completely between cycles. Therefore, we have chosen to include both IK,r and IK,s in our model, thus enabling it to simulate specific effects of antiarrhythmic drugs on IK,r, the buildup of IK,s at elevated heart rates, and the resulting changes, such as those in action potential duration (APD) and refractory period.
Our model of IK,r is based on data from Wang et al,44 Muraki et al,41 and Sanguinetti et al.45 Figs 6A⇓ and 6C⇓ show the steady-state activation (p̄a) and inactivation (p̄i) curves and the time constant of activation (τpa) compared with the available data. Since the inactivation of IK,r is very rapid compared with the activation, inactivation is modeled as being instantaneous. The expressions for IK,s are based on data recorded in human atrial myocytes by Wang and colleagues.42 44 Figs 6B⇓ and 6D⇓ show the characteristics of this current.
In order to verify our model of IK,r we chose to simulate the “ramp clamp” experiment of Muraki et al41 in rabbit atrial cells (their Fig 7⇓), in which the cell is subjected to an 0.8-V/s repolarizing ramp from an action potential peak potential down to −80 mV. This waveform is an approximation to the repolarization phase of an atrial action potential. The result (simulated response) is shown in Fig 7⇓, along with data from Muraki et al. Note that the simulated waveform is very similar, although in order to fit action potential data, we have had to reduce the size of IK,r to 15% of that shown in Fig 7⇓. Given that it is often difficult to detect any delayed rectifier current at all in most human atrial myocytes,2 it seems reasonable that this current should be assigned a very low density.
Time-Independent Currents: IK1, IB,Na, IB,Ca, INaK, INaCa, and ICaP
In the absence of reliable published data from human atrial cells for these currents, we have used expressions from the LMCG model13 with only minor scaling adjustments. One exception is the inward rectifier current, where we have found it necessary to make minor modifications to the rectifying characteristics in order to fit action potential data (slightly narrower outward “hump”). We have also adjusted the [K+]c dependence of IK1 to agree with data from our laboratory. Table 3⇓ lists the changes made in the time-independent currents compared with those in the LMCG model.
As in the LMCG model,13 we have included a fluid compartment formulation to monitor and account for changes in ion concentrations. These concentration changes can be a result of current flow across the cell membrane or of redistribution of ions within the cell (eg, uptake of Ca2+ by the SR or binding of Ca2+ to an intracellular buffer). Our fluid compartment model is similar to the one in the LMCG model. It includes descriptions of extracellular and intracellular spaces, formulations for Ca2+ uptake and release and the buffering action of calsequestrin, and troponin and calmodulin buffers in the intracellular medium. Compartment volumes and other ultrastructural properties, as well as expressions describing the binding of Ca2+ to intracellular troponin and calmodulin buffers and to calsequestrin in the SR release compartment, are identical to those of the LMCG model, except as noted in the following sections.
We have included a “cleft space” in our fluid compartment formulation, ie, a small restricted space surrounding the cell, in which accumulation or depletion of ions may occur (see Demir et al12 ). The cleft space is modeled as an unstirred fluid layer; ie, ions can be exchanged between the cleft space and the extracellular medium (in which all concentrations are assumed constant) only through diffusion as a result of a concentration gradient. Ratios between diffusion time constants for the ions involved (τNa, τK, and τCa) were calculated from values for ionic conductivity46 and the composition of the extracellular solution ([Cl−]o=140 mmol/L, [Na+]o=130 mmol/L, [K+]o=5.4 mmol/L, and [Ca2+]o=1.8 mmol/L). We have adjusted the size and diffusion properties of the cleft space so as to produce oscillations in cleft space [K+] ([K+]c) similar to experimental data.47 48
Electroneutral Na+ Influx
In order to achieve long-term stability in the ionic concentrations in the model, we have found it necessary to add a small (1.68-pA) electroneutral inward flux of Na+, denoted ΦNa,en. This flux could, for example, be accounted for in terms of electroneutral coupled transport mechanisms, such as Na+-H+ exchange and Na+-K+-2Cl cotransport. Modeling of these mechanisms, however, is beyond the scope of this work.
There are two major reasons for including this flux: First, the fact that long-term ionic homeostasis can be achieved with the addition of this small flux demonstrates that the sizes and other characteristics of the model elements are such that ionic homeostasis can reasonably be maintained. Second, if the ionic concentrations were allowed to change slowly from cycle to cycle (which would be the result if this flux were not included), the model would only be valid for short simulation times (seconds), for which this drift can safely be ignored. By ensuring long-term stability of the ionic concentrations, longer simulation times (minutes) become feasible and meaningful. Only with stable ionic concentrations can the model be used to simulate concentration changes as a result of rate changes or other interventions.
Our formulation for the SR is very similar to that of the LMCG model. However, we have made one important modification in accordance with recent evidence demonstrating that Ca2+ can accumulate in a small domain between the sarcolemma and the peripheral junctional SR and trigger Ca2+ release.27 49 Specifically, we have removed the voltage-dependent term in the formulation for activation of SR Ca2+ release and replaced it with a term dependent on Ca2+ concentration in the restricted subsarcolemmal domain, [Ca2+]d. The sole mechanism for SR Ca2+ release in our model is therefore Ca2+-induced Ca2+ release (CICR). Fig 1B⇑ includes a schematic representation of the model of the SR.
Stern50 has shown that in order for a CICR model to be stable, ie, capable of producing a response that is graded by the amount of Ca2+ that enters the cell through ICa,L, the trigger Ca2+ has to be separated from that released from the SR. Anatomically, this can be understood in terms of the concept of “release units” discussed, for example, by Isenberg and Han.51 According to this concept, Ca2+ release from the SR is recruited stepwise by the all-or-none activation of individual release units, consisting of one or more L-type Ca2+ channel and associated SR release channels. The activation of a release unit results in CICR within that unit only. The released Ca2+ then diffuses into the myoplasm, without directly affecting other units. This phenomenon has been incorporated in our model as a lumped mechanism, where the SR release channel senses [Ca2+] in the restricted domain ([Ca2+]d) but releases Ca2+ directly to the cytosol (Fig 1B⇑), thus separating trigger Ca2+ from that released from the SR.
A model of this type contains a large number of parameters that must be assigned values based on the available data. We have approached this part of the model development process in a two-step fashion, where the majority of the parameter values have been assigned in the first step, based on experimental studies of individual model components. This has the advantage that most of the parameters associated with an individual membrane current can be justified and assigned independently. Once the descriptions of individual membrane currents has been completed, one is left with a limited number of free parameters, most of which are scaling factors, such as the maximum conductance values for each ionic current. These remaining free parameters can then be determined using data for whole-cell responses (eg, action potential waveforms) or other constraints as indicated previously (eg, ionic homeostasis). It should be emphasized, however, that model development is very much an iterative process and that it has been necessary in some cases to modify individual current expressions to obtain acceptable fits to action potential data or (as in the case of INa) to use information from action potential recordings to resolve ambiguities in ion channel current data. The following sections describe the constraints and criteria used in order to assign values for the remaining free parameters.
The Quiescent Human Atrial Myocyte
In the absence of an external stimulus, a healthy human atrial cell is quiescent (does not contract or produce an action potential). In this quiescent state, the membrane potential comes to an equilibrium, “resting,” potential at which the net ion flux across the sarcolemma is zero. The resting membrane potential varies somewhat among individual cells, ranging from −70 to −80 mV. The resting membrane potential is the result of a precise balance between the time-independent background, pump, and exchanger currents (IK1, IB,Na, IB,Ca, INaK, ICaP, and INaCa). Although the resting state of the cell may seem less interesting than the active state during an action potential, an accurate description of the resting conditions is, in our experience, essential for successful modeling of the action potential. Moreover, the model of the resting state of the cell determines very important threshold characteristics and subthreshold properties, such as the input resistance of the cell. Very accurate simulation of these passive characteristics is essential before the cell model is to be used in distributed simulations, eg, of the propagation of electrical activity from one cell to another. We have modeled the resting state of the cell by adjusting the magnitudes of these currents so as to produce zero net transmembrane current at a resting potential of ≈75 mV and the following intracellular ion concentrations: [Na+]i≈8.5 mmol/L; [K+]i≈130.0 mmol/L, and [Ca2+]i≈60.0 nmol/L. The input resistance of the cell in this quiescent state is ≈600 MΩ, which agrees with our experimentally observed values.
The Activated Human Atrial Myocyte
When the cell is stimulated, it produces an action potential, the shape of which depends on the relative sizes of the ionic currents involved. Action potential data recorded from isolated human atrial myocytes are (in our experience) quite variable from cell to cell. It is therefore questionable whether it is possible to define the “normal human atrial action potential” in a meaningful way. We have chosen to fit our model to an action potential waveform that is representative of what is most often recorded from isolated human atrial cells in our laboratory. By establishing this “nominal model,” we have obtained a starting point from which the sensitivity of the action potential waveform to parameter perturbations may be studied (see “Results”). In addition to dictating the action potential waveform, the sizes (maximum conductance parameters) of the ionic currents also affect ionic homeostasis. The membrane currents involved in shaping the action potential therefore have to act in concert with those involved in the resting state to maintain constant ion concentrations at nominal stimulation rates. We have “tuned” our model so that ion concentrations remain constant from cycle to cycle at a stimulus frequency of 1 Hz.
There are two major aims of the simulations presented in the following sections: (1) to establish the validity and usefulness of our model by demonstrating that the expressions that are based on fits to voltage-clamp data for individual ionic currents also are able to accurately reconstruct action potential data and (2) to use the model to investigate the functional roles of different ionic currents, to study the sensitivity of the action potential waveform to the sizes of those currents, and, finally, to predict the whole-cell response to, for example, selected channel blocking drugs.
Simulated Action Potential Waveform
Fig 8A⇓ shows a simulated action potential waveform (solid line) compared with an action potential (dotted line) recorded at a temperature of 33°C and a stimulus frequency of 0.5 Hz. There is close agreement between the waveforms (the discrepancy at the beginning of the upstroke is due to a stimulation artifact that is not simulated). As mentioned previously, the action potential waveform varies significantly among individual cells and in multicellular preparations from the human atrium, presumably because of variations in the magnitudes of the underlying ionic currents. Nevertheless, the ability of this model to accurately reproduce a recorded action potential, in combination with the previously demonstrated fits to voltage-clamp data, lends credibility to the model.
In Fig 8⇑, panels B and C show the behavior of the membrane currents during an action potential. The first current to respond to a depolarizing stimulus pulse (delivered at time=100 ms) is INa, which activates rapidly, resulting in a very large but transient inward current. Note that INa is too large to be shown on the scale of Fig 8B⇑; its peak magnitude is ≈−5.8 nA, which corresponds to a maximum upstroke velocity of 116 V/s. INa is primarily responsible for the upstroke (phase 0) of the action potential, but as seen from Fig 8B⇑, a substantial amount of INa remains during the early peak phase of the action potential as a result of the second slower component of INa inactivation.
On depolarization of the cell, It, Isus, and ICa,L are also activated. However, It and Isus reach their peak magnitude faster than ICa,L, and their combined magnitude is thus larger than that of ICa,L early in the action potential. (This is because the peak of the action potential is close to the reversal potential for ICa,L.) The initial result on the action potential waveform is therefore a period of relatively rapid repolarization (phase 1), dominated by It. Since the time course of inactivation of ICa,L is slow compared with that of It, the net current gradually becomes dominated by ICa,L and Isus. A situation where the repolarizing effect of Isus (and the remaining It) is balanced by the depolarizing effect of ICa,L results. In the action potential waveform, the initial rapid repolarization (phase 1) is followed by a period during which the membrane potential levels off, forming a plateau (phase 2). Finally, as ICa,L slowly inactivates, the repolarizing effects of Isus become dominant, and the action potential enters its final repolarization phase (phase 3). During this phase, Isus is aided by the inward rectifier current, IK1, and the delayed rectifier currents, IK,r and IK,s, in repolarizing the cell membrane back to the resting potential (phase 4).
Simulated Ionic Fluxes
The fluid compartment part of this model monitors ion concentrations in the intracellular and cleft spaces. Valid modeling of the action potential requires not only the reconstruction of the action potential waveform but also a demonstration that this can be accomplished under conditions of ionic homeostasis at nominal heart rates. Table 4⇓ shows how our model has been tuned to achieve homeostasis at 1 Hz. The average charge transported across the sarcolemma for each ionic current has been computed by integrating each current over one cycle (1 s in the case of a quiescent, nonstimulated cell). Note that the sums of these average charges are zero for all ionic species at a stimulus rate of 1 Hz (Table 4⇓). When the cell is quiescent, there is a small net loss of intracellular Na+ and gain of intracellular K+. The existence of such an ionic imbalance at quiescence is supported by the observation by Bénardeau et al52 that trains of depolarizing pulses that activate INa can be used to hyperpolarize the resting potential of human atrial cells after a period of quiescence. As suggested by these authors, the hyperpolarization and stabilization of the resting potential may be caused by activation of the Na+-K+ pump after Na+ entry during the train of pulses. In our model, at a 2-Hz stimulus rate there will be a net gain of intracellular Na+ and loss of K+.
Fig 9⇓ illustrates the Ca2+ handling in the fluid compartment part of the model during the simulated action potential. A transient increase in [Ca2+]i occurs early in the action potential, raising [Ca2+]i from the very low diastolic levels (≈65 nmol/L) to a peak of ≈1.3 μmol/L (Fig 9B⇓). This rise in [Ca2+]i is primarily due to the rapid release of large amounts of Ca2+ from the SR (see Fig 9D⇓). Several processes are responsible for the decline of the intracellular Ca2+ transient. The most potent of these is the rapid binding of Ca2+ to the intracellular Ca2+ buffers (troponin and calmodulin), which is particularly important in shaping the early portions of the Ca2+ transient. As seen in Fig 9C⇓, the occupancies on these buffers increase rapidly as Ca2+ is released from the SR, thus “removing” large amounts of free Ca2+ from the cytosol. The uptake of Ca2+ by the SR also has a pronounced effect on the shape of the Ca2+ transient. This is also the primary pathway for actual removal of Ca2+ from the cytosol (as opposed to the “temporary storage” provided by the buffers), taking up intracellular Ca2+ as it dissociates from the buffers. In addition, some Ca2+ is removed from the cytosol via the Na+-Ca2+ exchanger, INaCa, and the Ca2+ pump, ICaP. As seen in Table 4⇑, INaCa and ICaP, on average, remove the amounts of Ca2+ that were brought into the cell via ICa,L and IB,Ca and thereby prevent a progressive buildup of cytosolic Ca2+ during repetitive stimulation.
Increasing the stimulus rate from the baseline (1 Hz) results in a change in ion concentrations in the intracellular medium (Table 4⇑) as well as in the cleft space surrounding the cell. For the ionic species that exist in relatively low concentrations in the extracellular medium (Ca2+ and K+), these changes can be significant. Fig 10⇓ shows how the intracellular and cleft space concentrations change when the stimulus rate is increased abruptly from 1 to 2 Hz. Note the progressive shift in [K+]c of ≈1 mmol/L.47 48 In contrast, the change in [Na+]c is negligible because of its high baseline value of 130 mmol/L. If the simulation in Fig 10⇓ is continued beyond the 20 s shown, [K+]c will reach a peak value of ≈6.3 mmol/L, after which it will begin to decline as a result of increased INaK activity due to increased [Na+]i and [K+]c (simulation not shown). The asymptotic values for [K+]c and [Na+]i are 5.6 and 10.0 mmol/L, respectively (reached after ≈10 minutes). This behavior is consistent with experimental observations.47
Parameter Sensitivity of the Action Potential Shape
As mentioned in the previous section, there is considerable variation in action potential shape among individual cells from the human atrium. Our working hypothesis is that many of these differences in action potential waveshape can be explained in terms of differences in the magnitudes of the ionic currents (caused by previous drug treatment and/or natural variability). In order to investigate possible mechanisms of action potential shape variability, it is therefore necessary to have an understanding of how changes in the magnitudes of different ionic currents affect the action potential shape. Such an understanding is equally important, of course, in identifying suitable “targets” for drug action aimed at modifying the action potential shape.
A valid mathematical model provides a method for this sensitivity analysis53 ; ie, it provides a method for studying how sensitive the state variables (eg, membrane voltage) are to perturbations in model parameters. We will restrict this analysis to a study of the sensitivity of the action potential waveform to changes in the sizes (maximum conductances) of the currents involved in shaping the action potential, although this type of analysis in principle can be used to study the sensitivity of any state variable to perturbations in any model parameter. Briefly, sensitivity analysis involves the computation of the partial derivative of the state variable of interest (in this case, membrane voltage) with respect to a model parameter. We have chosen to normalize these partial derivatives with respect to the nominal value of each parameter. Thus, we will present the results in terms of sensitivity functions, defined as: which can be interpreted as a proportionality factor relating a relative change in a model parameter to a resulting change in membrane voltage. The computational aspects of this method are outlined in Appendix B (see also Paulsen et al53 ).
Fig 11⇓ shows the result of the sensitivity analysis, ie, the sensitivity functions for the membrane voltage (action potential) with respect to the parameters of interest. Sensitivity functions were computed for all maximum conductance parameters in the model as well as for the scaling parameters for INaK and INaCa. The sensitivity functions for INa and IK,s have been omitted from Fig 11⇓, since these currents were found to have negligible influence on the action potential shape (small sensitivity functions), except for the obvious importance of INa during the upstroke. Several of the sensitivity functions have maximum absolute values of ≈50 mV during the late repolarization phase of the action potential. In other words, a 10% change in either one of these parameters would alter the membrane voltage by ≈5 mV during this phase of the action potential. Although this estimate is based on a linearization around the nominal parameter value and thus is most accurate for small perturbations, it can provide an indication of the approximate change expected for larger perturbations. Overall, we can anticipate that changes in these parameters in the ±50% range will produce significant changes in the action potential waveform.
Perhaps more important than the absolute sensitivity values, however, is the information provided by the time course of the sensitivity functions. As the ionic conductances change during the action potential, the sensitivity functions indicate which currents have the greatest influence on the action potential shape at each point in time. For example, it is clear that the action potential waveform is sensitive to ḡt primarily early in the peak phase of the action potential but that ḡsus and ḡCa,L rapidly become more important early in the plateau phase. Throughout the plateau phase, |εsus| (absolute value) and |εCa,L| are larger than |εt|, and this portion of the action potential waveform is therefore particularly sensitive to perturbations in ḡsus and ḡCa,L. Toward the end of repolarization, the action potential becomes very sensitive to the size of the time-independent currents involved in maintaining the resting potential. This analysis shows that (under the hypothesis that action potential changes can be explained in terms of changes in the magnitudes of ionic currents) ḡt, ḡsus, and ḡCa,L are the most important model parameters in determining the action potential waveshape during the peak and plateau phases of the action potential.
Roles of It and Isus in Repolarization
The effects on the human atrial action potential of blocking It with agents such as 4-aminopyridine, flecainide, and quinidine have been described in the literature.2 54 55 56 Partial block of It results in a slowing of the rate of repolarization of the action potential, particularly during the early repolarization phase (phase 1). This is, of course, consistent with the characteristics of It, and its role in the generation of the action potential as indicated by sensitivity analysis. In a recent study of the effects of some antiarrhythmic agents on It and Isus, Wang et al56 found that quinidine, in addition to blocking It, has a pronounced effect on Isus at clinically relevant concentrations. Considering that our sensitivity analysis indicates a prominent role for Isus in repolarization, the APD is expected to be quite sensitive to modulation of Isus magnitude. A detailed account of the role of It and Isus in the repolarization of the human atrial action potential is therefore essential for understanding the antiarrhythmic actions of quinidine.
Since pharmacological blocking agents used in experimental work usually affect more than one current and since their effects are often rate dependent, it is difficult to gain a quantitative understanding of the importance of a particular current (It or Isus) in repolarization from experimental results alone. In a computer model, however, it is possible to alter the characteristics of one ionic current in a controlled fashion, while leaving all other currents unaffected. Such simulations can be a valuable complement to experimental work. Given that drugs that prolong the APD (class III drugs) have been shown to be effective in the treatment of atrial arrhythmias,57 a thorough understanding of the influence of different ionic currents on the APD is needed.
Fig 12A⇓ shows the effects on the action potential of various degrees (30%, 60%, and 90%) of block of It. (In the present study, an x% block of It is simulated as an x% reduction in the maximum conductance, ḡt.) As observed experimentally, It block results in a broadening of the action potential peak during phase 1 of repolarization (refer to Fig 2 in Firek and Giles2 ). In addition, because of the elevation of the action potential peak and plateau levels, the contribution of Isus to repolarization is increased. The resulting prolongation of the APD is therefore only moderate, even for substantial reductions of It size. Fig 12B⇓ shows the effects of various degrees (15%, 30%, and 45%) of block of Isus on the action potential. In contrast to It block, inhibition of Isus primarily affects the plateau phase of the action potential, with little or no effect on the action potential peak. As a result, Isus block produces a more pronounced prolongation of the APD than does It block. For example, 30% block of Isus results in a 15% increase in APD at 90% repolarization (APD90) compared with the 5% increase in APD90 resulting from 30% It block. Many antiarrhythmic agents have effects on several ion channels. For example, according to Wang et al,56 quinidine blocks both It and Isus (in addition to its effects on Na+ channels). It is therefore of interest to study the effects on the action potential of combined It and Isus block. Fig 12C⇓ shows the result of a simulation in which both It and Isus have been reduced by 40% (approximately corresponding to the effect of 5 μmol/L quinidine at a stimulus rate of 1 Hz56). As expected, the result is essentially a combination of the previously demonstrated effects of It and Isus block, ie, a widened peak, an elevated action potential plateau, and a prolongation of APD90 of 27%.
Modulation of the Role of Isus by Baseline ICa,L, Isus, and IK,r Sizes
As discussed in “Model Development,” published data regarding the size of ICa,L and Isus are quite variable. Both these currents (as well as It) are known to be depressed in diseased human atrial cells22 26 and modulated by adrenergic stimulation.26 58 Furthermore, a recent study59 shows that the size of Isus is significantly reduced in cells obtained from patients in chronic atrial fibrillation compared with patients in normal sinus rhythm. It is therefore likely that a range of sizes of these two currents contributes to the physiological (and pathophysiological) behavior of the human atrial cell. Similarly, IK,r in our nominal model is very small, which is consistent with observations from our laboratory. Since results in other species60 indicate that this may be a consequence of the cell isolation techniques used for human atrial myocytes,2 this current may be considerably larger in vivo. Given these uncertainties in the actual sizes of several of the ionic currents, it is appropriate to investigate how the conclusions reached above are affected by our assumptions for these current sizes. We have chosen to focus on the role of Isus in repolarization, since this current is an important determinant of APD. Starting from our nominal model, we have performed a large number of simulations for different combinations of increased/decreased ICa,L, Isus, and IK,r. All simulations were performed at a stimulation rate of 1 Hz, and 20 cycles were allowed after each change of parameter values in order for any initial transient behavior to die out before the APD prolongation was evaluated. Fig 13⇓ shows how the APD prolongation resulting from a 50% reduction of Isus depends on the baseline sizes of ICa,L and IK,r. ICa,L and IK,r sizes are expressed as percentages of those in the nominal model; ie, the nominal model corresponds to 100% of both currents. A 7-fold (700%) increase in IK,r corresponds approximately to the size of IK,r observed in rabbit atrial myocytes.41 It is clear from Fig 13⇓ that the role of Isus as a major determinant of APD depends strongly on the size of IK,r current. The action potentials shown in the insets in Fig 13⇓ provide an indication of the underlying mechanism. When Isus is partially blocked, the action potential plateau is depolarized, which in turn increases the amount of IK,r (and IK,s) that is activated. This effect, which counteracts the APD-prolonging effect of Isus block, becomes stronger as the size of IK,r is increased. Similarly, Fig 14⇓ shows how the APD prolongation as a result of 50% Isus block depends on the baseline sizes of ICa,L and Isus. Again, the APD-prolonging effect of Isus block is strongly dependent on the baseline current densities. Generally, the APD prolongation as a result of partial Isus block becomes larger as the baseline size of ICa,L increases, provided that the increase in ICa,L is balanced by a comparable increase in Isus. Notice, however, in Fig 13⇓ as well as in Fig 14⇓, that when baseline ICa,L is increased without such a comparable increase in baseline outward current (Isus or IK,r), the APD-prolonging effect of Isus block is instead reduced. The underlying mechanism in this case is that the action potential plateau, both before and after Isus block, is prolonged and depolarized somewhat by the larger ICa,L. As a result, more IK,r (and IK,s) is activated, and the final repolarization phase becomes steeper and therefore less dependent on Isus.
Mathematical models form an important complement to experimental work in attempts to elucidate the ionic mechanisms underlying the action potential and other electrophysiological phenomena in cardiac cells. At a minimum, these models can provide a means of integrating data obtained from many different experiments and laboratories so that biophysically based explanations of complex nonlinear phenomena such as action potential initiation (excitation) and repolarization can be given. Mathematical models also provide a way of reviewing data in the context of the normal behavior of a cell (during an action potential), even when some of the data may have been obtained under completely different experimental conditions, eg, a voltage-clamp experiment. Moreover, as illustrated in the last few sections of “Results,” mathematical models can also have considerable predictive capabilities. After a model has been developed and carefully validated, it can be used to predict the response of the cell to selected drugs, experimental protocols, etc. Ideally, experimental work and model development should be carried out in close association, using the model to design and evaluate experiments and using experimental results to improve the model.
We have developed a mathematical model of the human atrial cell based primarily on data recorded from enzymatically isolated single human atrial cells. Our model is capable of accurately reconstructing a recorded human atrial action potential and illustrates the functional roles of the ionic currents. In addition, our model maintains ionic homeostasis at a nominal stimulus rate, demonstrating that the reconstruction of the action potential is accomplished using plausible current densities. We used the LMCG rabbit atrial model13 as a “starting point” for the model development. As a result, these two formulations are very similar in some aspects, particularly for the membrane currents where incomplete (or no) data from human atrial cells are available. However, there are some important differences between the electrophysiological responses of rabbit and human atrial cells; these provided the motivation for the development of this human atrial cell model. Perhaps most striking is the small rate dependence of It in human atrial cells compared with the very prominent rate dependence seen in rabbit atrial cells.13 38 This is mainly due to differences in the rates of recovery from inactivation of It in the two species. As discussed in “Membrane Currents,” It in human atrial cells recovers from inactivation much more rapidly than does It in rabbit atrial cells. Another important difference between human and rabbit atrial cells is the sustained outward current (denoted Isus in this article). Whereas this current in the rabbit atrium is believed to be carried mainly by Cl− ions,13 there is convincing evidence that Isus in human atrial cells is carried mainly by K+ ions.2 35 Both these differences are very important in studying the mechanisms of repolarization of the human atrial action potential and the effects of action potential–prolonging (class III) drugs that affect these currents.
Roles of It and Isus in Repolarization
Both It and Isus have important roles in the repolarization of the action potential of human atrium. In particular, Isus, because of its noninactivating characteristics, is necessary for repolarization to the resting potential. This important role in repolarization makes It and Isus potential targets for class III antiarrhythmic drugs, which are designed to prolong the APD. Indeed, a recent study of flecainide and quinidine,56 both known to prolong the human atrial action potential,54 55 shows that both these drugs produce a partial block of It. Quinidine also blocks Isus, which could explain its greater efficacy (compared with flecainide) in prolonging the APD.55 Our simulations of partial It and Isus block (see “Results”) produce prolongations of the APD that are comparable to experimentally observed effects of flecainide and quinidine, ie, a 27% prolongation of APD90 when It and Isus are both reduced by 40%. For comparison, Wang et al55 reported that 2.25 μmol/L quinidine increased APD95 by 33% in human atrial cells at a stimulation rate of 1 Hz. It should be noted that only some of the known effects of quinidine have been modeled; therefore, our results cannot be directly compared with these experimental observations. For example, the effects of quinidine on INa, as well as the state dependence of It block by quinidine,61 would have to be included in a more comprehensive treatment of quinidine effects. Nevertheless, the experimental observations agree very well with our model predictions and provide an independent “test” of how well our model describes the roles of It and Isus in repolarization.
Action potential generation involves a complex interaction among the ionic currents in a given cell type. The role of a particular ionic current in the action potential is therefore not determined solely by the characteristics of that current. We have investigated how the role of Isus in repolarization is affected when the baseline densities of ICa,L, Isus, and IK,r are varied within ranges that are relevant to the physiological and pathophysiological behavior of the human atrial cell. Our results demonstrate that Isus block will, in general, result in a prolongation of the action potential. The amount of prolongation, however, depends quite strongly on the baseline current densities. If the human atrial cell is assumed to have an IK,r density comparable to that observed in the rabbit atrium41 or if the ICa,L and Isus densities are reduced as observed in diseased cells,22 26 the APD prolongation resulting from Isus block may be considerably smaller than indicated by our nominal model. The efficacy of a drug targeting Isus would therefore be expected to depend critically on the disease state of the tissue. For example, based on the recent observation by Van Wagoner et al59 that Isus density is reduced in cells obtained from patients in chronic atrial fibrillation, the efficacy of an Isus-blocking drug may be limited in these patients.
Limitations of the Model
When using our model to gain insight into the electrophysiological responses of the human atrial cell, it is important to be aware of certain limitations, which are summarized by the following items:
1. The Hodgkin-Huxley formalism along with its concept of independent activation and inactivation “gating” variables has some important limitations. Notably, the processes of inactivation and recovery from inactivation (and analogously activation and deactivation) are governed by a single time constant. Experimental observations, however, often indicate that inactivation and recovery from inactivation occur with different time constants, even at the same membrane potential. In order to overcome this problem, one would have to use a more complicated modeling formalism that treats inactivation and recovery as two separate processes. To reduce the computational requirements of the model, we have chosen a “compromise” solution, in which time constant values are determined by measured inactivation kinetics at depolarized membrane potentials and by measured recovery kinetics at hyperpolarized potentials. In cases in which the time constants of inactivation and recovery are very different, this compromise results in unconventional time constant expressions, such as those for INa in Fig 2D⇑.
2. The available data regarding the intracellular Ca2+ transient and the Ca2+ handling in the SR have, with few exceptions,52 62 been recorded in cells from species other than humans. In addition, the understanding of the exact mechanisms involved in these phenomena is incomplete, and a quantitative model of SR Ca2+ release and uptake has not yet been developed. Our Ca2+-dependent formulation for SR Ca2+ release replaces the voltage-dependent formulation used in earlier models4 12 13 but is nevertheless only a qualitative description of this phenomenon. Other features of the SR function are the same as those in the LMCG model.13 Our SR formulation is therefore not based on human data.
3. The shape of the action potentials recorded from enzymatically isolated human atrial cells is variable, even at a fixed stimulus frequency. There are several reasons for this variability, including genuine heterogeneity between cells from different parts of the atrium, the disease states of the patients from whom the specimens are obtained, and the use of various drugs (eg, Ca2+ channel blockers and β-blockers) by the donors. It is important to acknowledge that this variability exists and that the exact action potential waveform in the model is chosen because it is representative of the action potential shapes that are most often observed in our single cell records. As indicated by sensitivity analysis, the shape of the human atrial action potential as described by our model is quite sensitive to variations in the strengths of three ionic currents (It, Isus, and ICa,L). Under the assumption that the currents involved and their kinetic properties are unchanged, our results therefore suggest that the experimentally observed variation in action potential shape is caused mainly by variations in these currents and can be explained within the framework of our model. However, it is also conceivable that regional diversity in the molecular basis of the currents in human atrium could give rise to regional differences in kinetics and pharmacological sensitivity. Most of the data on human atrial electrophysiology to date have been obtained from samples of the atrial appendage, and there are therefore little experimental data available regarding regional differences in human atrium. If and when data to suggest regional diversity become available, our model should provide a useful framework for predicting the consequences.
4. Our model of Ca2+-dependent inactivation of ICa,L uses a single lumped subsarcolemmal compartment in which Ca2+ accumulates and is therefore limited in its ability to simulate the effect of Ca2+ accumulation in restricted spaces close to each L-type Ca2+ channel. In contrast to voltage-dependent gating, where it is reasonable to assume that the controlling variable (membrane voltage) is spatially uniform (space-clamp conditions), this is in all likelihood not the case for subsarcolemmal Ca2+. Since inactivation is a nonlinear function of [Ca2+]d, it is not strictly correct to use a formulation in which all channels are subject to one (average) Ca2+ concentration.
Notwithstanding these limitations, this model provides the most complete description available of the ionic mechanisms underlying the human atrial action potential, and it is based on the available data. As a result, it provides a very useful tool for investigating fundamental electrophysiological responses of the human atrial cell, such as excitability, refractoriness, and the action of channel blocking drugs.
Tables 5 through 19⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓ contain all the equations, parameter values, and initial conditions necessary to carry out the simulations presented in this article. Unless otherwise noted, the units are as follows: time in seconds (s), voltage in millivolts (mV), concentration in millimoles/liter (mmol/L), current in picoamperes (pA), conductance in nanosiemens (nS), capacitance in nanofarads (nF), volume in nanoliters (nL), and temperature in kelvin (K). The stimulus used to evoke an action potential consists of a rectangular current pulse (Istim) with an amplitude of 280 pA and duration of 6 ms.
All simulations were performed by forward integration of the coupled system of differential equations using the CVODE solver package for ordinary differential equations. (CVODE was developed by S.D. Cohen and A.C. Hindmarsh at Lawrence Livermore National Laboratories, Livermore Calif) Sufficient accuracy was ensured by adjusting the temporal step size of the integration so that the local error in all state variables (as estimated by the CVODE algorithm) satisfied a relative error bound. Computer programs for the simulations were written in the C programming language under the UNIX operating system. Simulations were performed on Sun Microsystems Sparc workstations (Sparc 2, IPX) and on a Micron Millennia Pentium 166 PC running the Linux operating system. At a stimulus frequency of 1 Hz, one cycle (ie, 1 s of data) requires ≈0.9 s of CPU time on the Pentium 166 platform.
This Appendix contains a brief derivation of the method used to study the parameter sensitivity of the action potential waveform. The reader is referred to Reference 53 for more detailed information.
Starting from a system of the general form, we take the partial derivative of Equation 1 with respect to a parameter αk to obtain the following equation: The time derivative of the state vector x(α, t) can be written as follows: where the time derivative of α vanishes, since all parameters are assumed to be time invariant. Thus, from Equation 1, we have the following: Interchanging the order of integration, we obtain the following equation: where εk is the sensitivity of the state x to the parameter αk, ie, Note, finally, that Equation 5 defines a costate system53 of the same general form as Equation 1, ie, where The costate system for each parameter of interest can be integrated along with the original system, yielding the desired sensitivity functions.
τNa, τK, τCa INa
L-type Ca2+ current
Transient outward K+ current
Sustained outward K+ current
Slow delayed rectifier K+ current
Rapid delayed rectifier K+ current
Inwardly rectifying K+ current
Background Na+ current
Background Ca2+ current
Na+-K+ pump current
Sarcolemmal Ca2+ pump current
Na+-Ca2+ exchange current
Electroneutral Na+ influx
Ca2+ diffusion current from the diffusion-restricted subsarcolemmal space to the cytosol
Sarcoplasmic reticulum Ca2+ uptake current
Sarcoplasmic reticulum Ca2+ translocation current (from uptake to release compartment)
Sarcoplasmic reticulum Ca2+ release current
Na+ concentration in bulk (bathing) medium
K+ concentration in bulk (bathing) medium
Ca2+ concentration in bulk (bathing) medium
Na+ concentration in the extracellular cleft space
K+ concentration in the extracellular cleft space
Ca2+ concentration in the extracellular cleft space
Na+ concentration in the intracellular medium
K+ concentration in the intracellular medium
Ca2+ concentration in the intracellular medium
Mg2+ concentration in the intracellular medium
Ca2+ concentration in the restricted subsarcolemmal space
Ca2+ concentration in the sarcoplasmic reticulum uptake compartment
Ca2+ concentration in the sarcoplasmic reticulum release compartment
Equilibrium (Nernst) potential for Na+
Equilibrium (Nernst) potential for K+
Equilibrium (Nernst) potential for Ca2+
Apparent reversal potential for ICa,L (differs from ECa)
Permeability for INa
Maximum conductance for ICa,L
Maximum conductance for It
Maximum conductance for Isus
Maximum conductance for IK,s
Maximum conductance for IK,r
Maximum conductance for IK1
Maximum conductance for IB,Na
Maximum conductance for IB,Ca
Activation gating variable for INa
Fast and slow inactivation gating variables for INa
Activation gating variable for ICa,L
Fast and slow inactivation gating variables for ICa,L
[Ca2+]d-dependent ratio of fast (fL1) to slow (fL2) inactivation of ICa,L
Half-maximum Ca2+ binding concentration for fCa
Activation gating variable for It
Inactivation gating variable for It
Rapidly and slowly recovering inactivation gating variables for It
Activation gating variable for Isus
Inactivation gating variable for Isus
Activation gating variable for IK,s
Activation gating variable for IK,r
Inactivation gating variable (instantaneous) for IK,r
¯m, ¯h1, …
Steady-state value of m, h1, etc
Relative amount of “inactive precursor” in the Irel formulation
Relative amount of “activator” in the Irel formulation
Activation time constant for INa
Fast and slow inactivation time constants for INa
Activation time constant for ICa,L
Fast and slow inactivation time constants for ICa,L
Activation time constant for It
Inactivation time constant for It
Rapidly and slowly recovering inactivation time constants for It
Activation time constant for Isus
Inactivation time constant for Isus
Activation time constant for IK,s
Activation time constant for IK,r
Fractional occupancy of the calmodulin buffer by Ca2+
Fractional occupancy of the troponin-Ca2+ buffer by Ca2+
Fractional occupancy of the troponin-Mg2+ buffer by Ca2+
Fractional occupancy of the troponin-Mg2+ buffer by Mg2+
Fractional occupancy of the calsequestrin buffer (in the sarcoplasmic reticulum release compartment) by Ca2+
Universal gas constant
Volume of the extracellular cleft space
Total cytosolic volume
Volume of the diffusion-restricted subsarcolemmal space
Volume of the sarcoplasmic reticulum uptake compartment
Volume of the sarcoplasmic reticulum release compartment
τNa, τK, τCa
Time constant of diffusion of Na+, K+, and Ca2+ from the bulk medium to the extracellular cleft space
Time constant of diffusion from the restricted subsarcolemmal space to the cytosol
Maximum Na+-K+ pump current
Half-maximum K+ binding concentration for INaK
Half-maximum Na+ binding concentration for INaK
Maximum Ca2+ pump current
Half-maximum Ca2+ binding concentration for ICaP
Scaling factor for INaCa
Position of energy barrier controlling voltage dependence of INaCa
Denominator constant for INaCa
Maximum sarcoplasmic reticulum uptake current
Half-maximum binding concentration for [Ca2+]i to Iup
Half-maximum binding concentration for [Ca2+]up to Iup
Ratio of forward to back reactions for Iup
Time constant of diffusion (“translocation”) of Ca2+ from sarcoplasmic reticulum uptake to release compartment
Scaling factor for Irel
Half-activation [Ca2+]i for Irel
Half-activation [Ca2+]d for Irel
Recovery rate constant for the sarcoplasmic reticulum release channel
Mr Nygren was supported in part by a stipend from a Canadian Medical Research Council group grant. Drs Fiset and Firek are grateful for the support of Alberta Heritage Foundation for Medical Research fellowships. Dr Giles is supported as a medical scientist by the Alberta Heritage Foundation for Medical Research and receives ongoing support from the Canadian Medical Research Council and the Heart and Stroke Foundation of Canada.
- Received June 9, 1997.
- Accepted September 30, 1997.
- © 1998 American Heart Association, Inc.
Shibata EF, Drury T, Refsum H, Aldrete V, Giles W. Contributions of a transient outward current to repolarization in human atrium. Am J Physiol. 1989;257:H1773–H1781.
Firek L, Giles WR. Outward currents underlying repolarization in human atrial myocytes. Cardiovasc Res. 1995;39:31–38.
DiFrancesco D, Noble D. A model of cardiac electrical activity incorporating ionic pumps and concentration changes. Philos Trans R Soc Lond B Biol Sci. 1985;307:353–398.
Rasmusson RL, Clark JW, Giles WR, Robinson K, Clark RB, Shibata EF, Campbell DL. A mathematical model of electrophysiological activity in a bullfrog atrial cell. Am J Physiol. 1990;259:H370–H389.
Rasmusson RL, Clark JW, Giles WR, Shibata EF, Campbell DL. A mathematical model of a bullfrog cardiac pacemaker cell. Am J Physiol. 1990;259:H352–H369.
Luo C-H, Rudy Y. A model of the ventricular cardiac action potential: depolarization, repolarization and their interaction. Circ Res. 1991;68:1501–1526.
Luo C-H, Rudy Y. A dynamic model of the cardiac ventricular action potential, I: simulations of ionic currents and concentration changes. Circ Res. 1994;74:1071–1096.
Luo C-H, Rudy Y. A dynamic model of the cardiac ventricular action potential, II: afterdepolarizations, triggered activity, and potentiation. Circ Res. 1994;74:1097–1113.
Zeng J, Laurita KR, Rosenbaum DS, Rudy Y. Two components of the delayed rectifier K+ current in ventricular myocytes of the guinea pig type: theoretical formulation and their role in repolarization. Circ Res. 1995;77:140–152.
Demir SS, Clark JW, Murphey CR, Giles WR. A mathematical model of a rabbit sinoatrial node cell. Am J Physiol. 1994;266:C832–C852.
Lindblad DS, Murphey CR, Clark JW, Giles WR. A model of the action potential and underlying membrane currents in a rabbit atrial cell. Am J Physiol. 1996;271:H1666–H1691.
Ilnicki T. Electrophysiological and Mechanical Measurements in Human and Rabbit Atria. Calgary, Canada: University of Calgary; 1987. Thesis.
Sakakibara Y, Wasserstrom JA, Furukawa T, Jia H, Arentzen CE, Hartz RS, Singer DH. Characterization of the sodium current in single human atrial myocytes. Circ Res. 1992;71:535–546.
Sakakibara Y, Furukawa T, Singer, Donald H, Jia H, Backer CL, Arentzen CE, Wasserstrom JA. Sodium current in isolated human ventricular myocytes. Am J Physiol. 1993;265:H1301–H1309.
Wendt DJ, Starmer F, Grant AO. Na channel kinetics remain stable during perforated-patch recordings. Am J Physiol. 1992;263:C1234–C1240.
Le Grand B, Hatem S, Deroubaix E, Couétil J-P, Coraboeuf E. Calcium current depression in isolated human atrial myocytes after cessation of chronic treatment with calcium antagonists. Circ Res. 1991;69:292–300.
Li GR, Nattel S. Properties of human atrial ICa at physiological temperatures and relevance to action potential. Am J Physiol. 1997;272:H227–H235.
Sun XH, Protasi F, Takahashi M, Takeshima H, Ferguson DG, Franzini-Armstrong C. Molecular architecture of membranes involved in excitation-contraction coupling of cardiac muscle. J Cell Biol. 1995;129:659–671.
Giles WR, Clark RB, Braun AP. Ca2+-independent transient outward current in mammalian heart. In: Morad M, Ebashi S, Trautwein W, Kurachi Y, eds. Molecular Physiology and Pharmacology of Cardiac Ion Channels and Transporters. Amsterdam, the Netherlands: Kluwer Publishing Ltd; 1996:141–168.
Escande D, Coulombe A, Faivre JF, Deroubaix E, Coraboeuf E. Two types of transient outward currents in adult human atrial cells. Am J Physiol. 1987;252:H142–H148.
Wettwer E, Amos G, Gath J, Zerkowski HR, Reidemeister JC, Ravens U. Transient outward current in human and rat ventricular myocytes. Cardiovasc Res. 1993;27:1662–1669.
Fedida D, Wible B, Wang Z, Fermini B, Faust F, Nattel S, Brown AM. Identity of a novel delayed rectifier current from human heart with a cloned K+ channel current. Circ Res. 1993;73:210–216.
Wang Z, Fermini B, Nattel S. Sustained depolarization-induced outward current in human atrial myocytes: evidence for a novel delayed rectifier K+ current similar to Kv1.5 cloned channel currents. Circ Res. 1993;73:1061–1076.
Näbauer M, Beuckelmann DJ, Erdmann E. Characteristics of transient outward current in human ventricular myocytes from patients with terminal heart failure. Circ Res. 1993;73:386–394.
Agus ZS, Dukes ID, Morad M. Divalent cations modulate the transient outward current in rat ventricular myocytes. Am J Physiol. 1991;261:C310–C318.
Fermini B, Wang Z, Duan D, Nattel S. Differences in rate dependence of transient outward current in rabbit and human atrium. Am J Physiol. 1992;263:H1747–H1754.
Muraki K, Imaizumi Y, Watanabe M, Habuchi Y, Giles WR. Delayed rectifier K+ current in rabbit atrial myocytes. Am J Physiol. 1995;269:H524–H532.
Wang Z, Fermini B, Nattel S. Delayed rectifier outward current and repolarization in human atrial myocytes. Circ Res. 1993;73:276–285.
Lide DR, ed. CRC Handbook of Chemistry and Physics. Cleveland, Ohio: CRC Press; 1992.
Kunze DL. Rate-dependent changes in extracellular potassium in the rabbit atrium. Circ Res. 1977;41:122–127.
Cohen I, Kline R. K+ fluctuations in the extracellular spaces of cardiac muscle: evidence from the voltage-clamp and extracellular K+-selective microelectrodes. Circ Res. 1982;50:1–16.
Bénardeau A, Hatem SN, Rucker-Martin C, Le Grand B, Mace L, Dervanian P, Mercadier JJ, Coraboeuf E. Contribution of Na+/Ca2+ exchange to action potential of human atrial myocytes. Am J Physiol. 1996;271:H1151–H1161.
Wang Z, Pelletier LC, Talajic M, Nattel S. Effects of flecainide and quinidine on human atrial action potentials: role of rate-dependence and comparison with guinea pig, rabbit, and dog tissues. Circulation. 1990;82:274–283.
Wang Z, Fermini B, Nattel S. Effects of flecainide, quinidine, and 4-aminopyridine on transient outward and ultrarapid delayed rectifier currents in human atrial myocytes. J Pharmacol Exp Ther. 1995;272:184–196.
Kottkamp H, Haverkamp W, Borggrefe M, Breithardt G. The role of class III antiarrhythmic drugs in atrial fibrillation. In: Olsson SB, Allessie MA, Campbell RWF, eds. Atrial Fibrillation: Mechanisms and Therapeutic Strategies. Armonk, NY: Futura Publishing Co Inc; 1994:287–306.
Li G-R, Feng J, Wang Z, Fermini B, Nattel S. Adrenergic modulation of ultrarapid delayed rectifier K+ current in human atrial myocytes. Circ Res. 1995;78:903–915.
Van Wagoner DR, Pond AL, McCarthy PM, Trimmer JS, Nerbonne JM. Outward K+ current densities and Kv1.5 expression are reduced in chronic human atrial fibrillation. Circ Res. 1997;80:772–781.
Yue L, Feng J, Li GR, Nattel S. Transient outward and delayed rectifier currents in canine atrium: properties and role of isolation methods. Am J Physiol. 1996;270:H2157–H2168.
Hatem SN, Bénardeau A, Rücker-Martin C, Marty I, de Chamisso P, Villaz M, Mercadier J-J. Different compartments of sarcoplasmic reticulum participate in the excitation-contraction coupling process in human atrial myocytes. Circ Res. 1997;80:345–353.