## Abstract

Dogs have been used extensively to study atrial arrhythmias, but there are no published mathematical models of the canine atrial action potential (AP). To obtain insights into the ionic mechanisms governing canine atrial AP properties, we incorporated formulations of K^{+}, Na^{+}, Ca^{2+}, and Cl^{−} currents, based on measurements in canine atrial myocytes, into a mathematical model of the AP. The rate-dependent behavior of model APs corresponded to experimental measurements and pointed to a central role for L-type Ca^{2+} current inactivation in rate adaptation. Incorporating previously described regional ionic current variations into the model largely reproduced AP forms characteristic of the corresponding right atrial regions (appendage, pectinate muscle, crista terminalis, and atrioventricular ring). When ionic alterations induced by tachycardia-dependent remodeling were incorporated, the model reproduced qualitatively the AP features constituting the cellular substrate for atrial fibrillation. We conclude that this ionic model of the canine atrial AP agrees well with experimental measurements and gives potential insights into mechanisms underlying functionally important electrophysiological phenomena in canine atrium.

- action potential duration
- atrial fibrillation
- ion channels
- rate adaptation
- regional heterogeneity
- mathematical model

Atrial arrhythmiassuch as atrial fibrillation (AF) represent common clinical problems that remain difficult to treat. Experimental dog models have been used widely to study atrial arrhythmia mechanisms in vivo (26). In vivo studies have provided useful insights into the role of atrial tachycardia-induced remodeling (36), electrical heterogeneity (10), and pathological alterations (3,4) in promoting the occurrence of atrial reentrant arrhythmias. More recently, patch clamp studies in isolated canine atrial myocytes have clarified the properties of a variety of ionic currents in dog atrium (38-40) and suggested ionic mechanisms for regional variations in action potential (AP) properties (11) and tachycardia-induced AP remodeling (38). It is of great interest and potential importance to link these two levels of observation to determine the ionic mechanisms that control the occurrence of arrhythmias in vivo.

Pharmacological tools can be used to evaluate the roles of individual ionic currents (38-40). However, the lack of specificity of pharmacological probes is an important limitation. In addition, changes in one ionic current alter the voltage-time trajectory of the AP and therefore secondarily affect other currents. Therefore, the changes in the AP caused by even a specific ion current inhibitor cannot be interpreted as reflecting the role of only that current. One approach to obtaining a more precise appreciation of the ionic determinants of the AP is to create a mathematical model incorporating realistic biophysical simulations of the ionic processes believed to be involved. The model can be interrogated to assess its agreement with physiological behavior in a variety of conditions to determine its applicability, and the role of individual components can be tested by modifying them and evaluating the resulting impact on the AP. The resulting simulations can then be compared to directly measured changes under experimental conditions.

The first mathematical model of the AP was developed by Hodgkin and Huxley (16) to simulate the electrical behavior of the squid giant axon. Since then, mathematical models of APs based on formulations of ionic currents, pumps, and exchangers have provided insights into properties of rabbit atrial (21) and sinoatrial node (8), guinea pig ventricular (22), bullfrog atrial (29), Purkinje fiber (24), and canine ventricular (37) APs. More recently, models of the human atrial AP have been published (6,27). These models can account for a wide range of important behaviors and have been used to analyze the effects on human atrial electrophysiology of tachycardia-induced remodeling and selective K^{+} channel blockade (7).

There is no published mathematical model of the canine atrial AP. Such a model would be valuable to interpret better the extensive information available about atrial arrhythmias in vivo in the dog and to consolidate the increasing body of knowledge regarding canine atrial ionic mechanisms. The objectives of the present study were*1*) to develop a mathematical model of the canine atrial AP using information obtained from the direct measurement of ionic currents in canine atrial myocytes, and *2*) to evaluate the agreement between model APs and APs from different experimental paradigms including varying activation rate, discrepant locations of recording in the atrium, and electrical remodeling caused by long-term atrial tachycardia.

*R*- Gas constant
*T*- Temperature
*F*- Faraday constant
*C*_{m}- Membrane capacitance
- AP
- Action potential
*V*_{max}- Maximal AP upstroke velocity
- APA
- AP amplitude
- APO
- AP overshoot
- APD
- AP duration
- APD
_{90} - APD to 90% repolarization
- APD
_{95} - APD to 95% repolarization
- AF
- Atrial fibrillation
- SR
- Sarcoplasmic reticulum
- JSR
- Junctional SR, SR release compartment
- NSR
- Network SR, SR uptake compartment
- V
_{cell} - Cell volume
- V
_{i} - Intracellular volume
- V
_{up} - NSR volume
- V
_{rel} - JSR volume
- [
*X*]_{o} - Extracellular concentration of ion
*X* - [
*X*]_{i} - Intracellular concentration of ion
*X* - Cmdn
- Calmodulin, sarcoplasmic Ca
^{2+}buffer - Trpn
- Troponin, sarcoplasmic Ca
^{2+}buffer - Csqn
- Calsequestrin, JSR Ca
^{2+}buffer - [Ca
^{2+}]_{rel} - Ca
^{2+}concentration in JSR - [Ca
^{2+}]_{up} - Ca
^{2+}concentration in NSR - [Ca
^{2+}]_{Cmdn} - Ca
^{2+}-bound calmodulin concentration - [Ca
^{2+}]_{Trpn} - Ca
^{2+}-bound troponin concentration - [Ca
^{2+}]_{Csqn} - Ca
^{2+}-bound calsequestrin concentration *E*_{X}- Equilibrium potential for ion
*X* *I*_{ion}- Total sarcolemmal ionic current
*I*_{stim}- Stimulus current
- α
_{x} - Forward rate constant for gating variable
*x* - β
_{x} - Reverse rate constant for gating variable
*x* - τ
_{x} - Time constant for gating variable
*x* *x*_{∞}- Steady-state relation for gating variable
*x* *I*_{Na}- Fast inward Na
^{+}current *g*_{Na}- Maximal
*I*_{Na}conductance *m**I*_{Na}activation variable*h**I*_{Na}fast inactivation variable*j**I*_{Na}slow inactivation variable*I*_{K1}- Inward rectifier K
^{+}current *g*_{K1}- Maximal
*I*_{K1}conductance *I*_{to}- Transient outward K
^{+}current *g*_{to}- Maximal
*I*_{to}conductance *o*_{a}*I*_{to}activation variable*o*_{i}*I*_{to}inactivation variable*I*_{Kur,d}- Ultrarapid delayed rectifier K
^{+}current *g*_{Kur,d}- Maximal
*I*_{Kur,d}conductance *u*_{a}*I*_{Kur,d}activation variable*u*_{i}*I*_{Kur,d}inactivation variable*I*_{Kr}- Rapid delayed rectifier K
^{+}current *g*_{Kr}- Maximal
*I*_{Kr}conductance *x*_{r}*I*_{Kr}activation variable*I*_{Ks}- Slow delayed rectifier K
^{+}current *g*_{Ks}- Maximal
*I*_{Ks}conductance *x*_{s}*I*_{Ks}activation variable*I*_{Ca}- Inward Ca
^{2+}current *g*_{Ca}- Maximal
*I*_{Ca}conductance *d**I*_{Ca}activation variable*f**I*_{Ca}voltage-dependent inactivation variable*f*_{Ca}*I*_{Ca}Ca^{2+}-dependent inactivation variable*I*_{NaCa}- Na
^{+}/Ca^{2+}exchanger current *I*_{NaCa(max)}*I*_{NaCa}scaling factor*I*_{p,Ca}- Sarcoplasmic Ca
^{2+}pump current *I*_{p,Ca(max)}- Maximal
*I*_{p,Ca} *I*_{NaK}- Na
^{+}-K^{+}pump current *f*_{NaK}- Voltage-dependence parameter for
*I*_{NaK} - ς
- [Na
^{+}]_{o}-dependence parameter for*I*_{NaK} *K*_{m,Na(i)}- [Na
^{+}]_{i}half-saturation constant for*I*_{NaK} *K*_{m,K(o)}- [K
^{+}]_{o}half-saturation constant for*I*_{NaK} *I*_{Cl,Ca}- Ca
^{2+}-activated Cl^{−}current *q*_{Ca}- Ca
^{2+}flux-dependent activation gating variable for*I*_{Cl,Ca} *K*_{m,Na}- [Na
^{+}]_{o}saturation constant for*I*_{NaCa} *K*_{m,Ca}- [Ca
^{2+}]_{o}saturation constant for*I*_{NaCa} *k*_{sat}- Saturation constant for
*I*_{NaCa} *I*_{b,Na}- Background Na
^{+}current *g*_{b,Na}- Maximal
*I*_{b,Na}conductance *I*_{b,Ca}- Background Ca
^{2+}current *g*_{b,Ca}- Maximal
*I*_{b,Ca}conductance *I*_{rel}- Ca
^{2+}release current from the JSR *k*_{rel}- Maximal Ca
^{2+}release rate for*I*_{rel} *u*- Activation gating variable for
*I*_{rel} *v*- Ca
^{2+}flux-dependent inactivation gating variable for*I*_{rel} *w*- Voltage-dependent inactivation gating variable for
*I*_{rel} - F
_{n} - Sarcoplasmic Ca
^{2+}flux signal for*I*_{rel} *I*_{up}- Ca
^{2+}uptake current into the NSR *I*_{up(max)}- Maximal Ca
^{2+}uptake rate for*I*_{up} - [Ca
^{2+}]_{up(max)} - Maximal Ca
^{2+}concentration in NSR *I*_{tr}- Ca
^{2+}transfer current from NSR to JSR - τ
_{tr} - Ca
^{2+}transfer time constant *I*_{up,leak}- Ca
^{2+}leak current from the NSR - [Cmdn]
_{max} - Total calmodulin concentration in myoplasm
- [Trpn]
_{max} - Total troponin concentration in myoplasm
- [Csqn]
_{max} - Total calsequestrin concentration in JSR
*z*- Valence

## MODEL DESCRIPTION

The cell membrane was modeled as a capacitor connected in parallel with variable resistances (ion channels) and batteries (driving forces) following the Hodgkin-Huxley formalism for an excitable membrane (16). In this formalism, each ionic current is proportional to the product of the driving force and the appropriate membrane conductance, which is in turn governed by gating functions with activation and inactivation variables. The activation and inactivation variables are represented as probabilities of ion channel opening, varying between 0 and 1, and are generally voltage and time dependent. The rate of change in the transmembrane potential (*V*) is given by
Equation 1where *I*
_{ion} and *I*
_{stim}are the total transmembrane ionic current and stimulus current, respectively, and *C*
_{m} is the total membrane capacitance. The total transmembrane ionic current is given by
Equation 2
where *I*
_{Na} is the fast Na^{+} current, *I*
_{K1} is the inward rectifier K^{+} current, *I*
_{to} is the transient outward K^{+} current, *I*
_{Kur,d}is the dog ultrarapid delayed rectifier K^{+} current,*I*
_{Kr} and *I*
_{Ks} are the rapid and slow delayed rectifier K^{+} current components, respectively, *I*
_{Ca} is the L-type Ca^{2+} current, *I*
_{Cl,Ca} is the Ca^{2+}-activated Cl^{−} current,*I*
_{NaCa} is the Na^{+}/Ca^{2+} exchanger current,*I*
_{NaK} is the Na^{+}-K^{+} pump current, and*I*
_{b,Na} and *I*
_{b,Ca} are the background Na^{+} and Ca^{2+} currents, respectively.

The model constantly monitors intracellular concentrations of Na^{+}, K^{+}, Ca^{2+}, and Cl^{−}and maintains constant extracellular concentrations. Figure1 is a schematic representation of the ionic components of the model, which include membrane currents, pumps, and exchangers. The intracellular space includes network (NSR) and junctional compartments of the sarcoplasmic reticulum (JSR) that play a role in the intracellular handling of Ca^{2+} (uptake from and release to the myoplasm, respectively). Unless otherwise noted, physical units are as follows: time (*t*) is in milliseconds,*V* is in millivolts, *C*
_{m} is in picofarads, current densities are in picoamperes per picofarad, conductances (*g*) are in nanosiemens per picofarad, and concentrations are in millimoles per liter.

The computer software encoding the model was written in C with the use of double-precision arithmetic. All simulations were performed on a Dell OptiPlex GX1 personal computer with an Intel Pentium II processor at 450 Hz using a RedHat Linux operating system. Time derivatives were integrated with a modified Euler method employing fixed time steps of 0.005 ms. This method ensured that the largest change in transmembrane potential over a single time iteration did not exceed 0.25 mV (33). To ensure convergence of solutions using this approach, simulations were also obtained with time steps of 0.001 ms (reduced by 80%). The differences between the APs obtained with time steps of 0.001 and 0.005 ms varied between 0.0031 ± 0.016 mV (mean ± SD) and did not exceed 0.23 mV at any point during the AP. The full system of equations is given in the , and model constants are given in Table 1.

### Membrane Currents

#### Fast Na^{+}* current.*

The model implemented *I*
_{Na} using the modification of the Ebihara-Johnson model (9) proposed by Luo and Rudy (22) and applied previously in our human atrial AP model (6). Model *I*
_{Na} is given by
Equation 3where *g*
_{Na} is the maximal Na^{+} conductance, *m* is the activation variable, and *h* and *j* are the inactivation variables.

#### Transient outward K^{+}* current.*

We formulated *I*
_{to} on the basis of results obtained from our laboratory in isolated canine atrial myocytes (38, 39). Figure2
*A* shows steady-state values for *I*
_{to} inactivation and activation in the model. Experimental values from Yue et al. (39) were fit to the functions given in the
by using a nonlinear least-squares algorithm (23). Figure2
*B* shows corresponding model and experimental time constants. The activation time constants were bell shaped, with a maximum of 5.9 ms at −27 mV. Experimental values for onset and removal of inactivation are represented by open and filled symbols, respectively. The current is given by
Equation 4where *g*
_{to} is the maximal*I*
_{to} conductance, *o*
_{a} is the activation variable, and *o*
_{i} is the inactivation variable. Figure 2
*C* shows the simulated response of *I*
_{to} to a voltage-clamp pulse protocol (*inset*). The *g*
_{to} was adjusted to obtain current amplitudes that agreed with experimental observations. The peak current-voltage (*I-V*) relationship (Fig. 2
*D*) obtained by the model shows excellent agreement with the experimental values reported by Yue et al. (39).

#### Dog ultrarapid delayed rectifier K^{+}*current.*

The current is given by
Equation 5
Equation 6where *g*
_{Kur,d} is the voltage-dependent maximal conductance, and *u*
_{a}and *u*
_{i} are the activation and inactivation variables, respectively. Figure 3,*A* and *B*, shows steady-state values and time constants for activation and inactivation of*I*
_{Kur,d}. The activation time constants shown in Fig. 3
*B* were fit to values reported by Yue et al. (40). Inactivation of *I*
_{Kur,d} is very slow. In the absence of canine-specific measurements for*I*
_{Kur,d} inactivation at 37°C, we used the inactivation steady-state values and time constants for human atrial myocytes reported by our laboratory (34) and used previously in our human AP model (6). Figure 3
*C*shows *I*
_{Kur,d} elicited by a simulated voltage-clamp pulse protocol. The peak *I-V* relationship in Fig. 3
*D* shows good agreement with data from Yue et al. (40).

#### Rapid and slow K^{+}* delayed rectifier current components.*

The currents are given by
Equation 7
Equation 8where *g*
_{Kr} and *g*
_{Ks}are the maximal current densities and *x*
_{r} and*x*
_{s} are the activation variables for the respective current components. Voltage-dependent activation steady states are shown in Fig. 4
*A*; experimental values are from Li et al. (19). Figure4
*B* shows the activation time constants along with experimental data. Figure 4, *C* and *D*, shows model currents elicited by the pulse protocol (*inset*). Figure 4,*E* and *F*, shows the peak *I-V*relationships provided by the model and obtained from the experimental data of Li et al. (19).

#### Inward rectifier K^{+}* current.*

The model follows the formulation of *I*
_{K1} in our human AP model (6). Maximal conductance (*g*
_{K1}) was adjusted to obtain a peak*I-V* relationship corresponding to experimental results (38). The voltage-dependent, time-independent current is given by
Equation 9

#### Ca^{2+}* current.*

We based our formulation of *I*
_{Ca} on experimental data (38) and previous AP models (6, 22). The current is given by
Equation 10where *g*
_{Ca} is the maximal conductance. The formulation includes voltage-dependent activation (*d*) and voltage- and Ca^{2+}-dependent inactivation (*f* and *f*
_{Ca}). The computed equilibrium potential for Ca^{2+} (*Eq. EA10
*) is between 100 and 130 mV, which is inconsistent with experimental data for *I*
_{Ca} reversal (11, 19, 38), in accordance with previous observations. The reversal potential for*I*
_{Ca} was therefore fixed at 65 mV as in previous models (6, 21). The steady-state activation and inactivation variables (Fig.5
*A*) and the time constants for voltage-dependent inactivation (τ_{f}) (Fig. 5
*B*) were fit to experimental values (38). In the absence of reported activation time constants specific to canine atria, we adopted the formulation of Luo and Rudy (22). Figure 5
*C*shows steady-state Ca^{2+}-dependent*I*
_{Ca} inactivation. Figure 5
*E* shows*I*
_{Ca} elicited by a simulated pulse protocol. The peak *I-V *relationship obtained from the model (Fig.5
*F*) is U-shaped, reversing at 65 mV, and agrees with experimental data (38). Ca^{2+} transients elicited at different rates are shown in Fig. 5
*D*.

#### Ca^{2+}*-activated Cl*^{−}current.

^{−}current.

The detailed mechanism of *I*
_{Cl,Ca} activation remains poorly understood, although direct activation by free cytoplasmic Ca^{2+} appears central. Factors influencing*I*
_{Cl,Ca} include channel density, colocalization with sarcolemmal Ca^{2+} channels, Ca^{2+} load in the sarcoplasmic reticulum (SR), and Ca^{2+}-induced Ca^{2+} release (5). In the model, activation of*I*
_{Cl,Ca} was initiated by Ca^{2+} flux into the cell, determining the increase in intracellular Ca^{2+} concentration ([Ca^{2+}]_{i}) in a localized region between the closely coupled sarcolemmal*I*
_{Ca} and the Ca^{2+} release channels of the SR (18, 41). The current is given by
Equation 11where *g*
_{Cl,Ca} is the maximal current conductance and *q*
_{Ca} represents the Ca^{2+} flux-dependent activation variable as defined in Fig.6
*A*. In response to depolarizing pulses, model-derived *I*
_{Cl,Ca} (Fig.6
*B*) decayed fully within 20 ms of the onset of depolarization. The peak *I-V* relationship (Fig.6
*C*) resembles experimental values reported by Yue et al. (38).

#### Na^{+}*-K*^{+}pump current and Na^{+}/Ca^{2+}exchanger current.

^{+}pump current and Na

^{+}/Ca

^{2+}exchanger current.

We implemented the Na^{+}-K^{+} pump (Eqs.*
EA57-EA59
*) by following previous models (6,22). The parameter for maximal current [*I*
_{NaK(max)}] was adjusted to maintain stable intracellular ion concentrations at rest. The exchanger current formulation (*Eq. EA60
*) similarly followed previously published formulations (6, 22).

#### Background Na^{+}* and Ca*^{2+} currents.

^{2+}currents.

The amplitudes of the background currents were adjusted to maintain stable resting intracellular ion concentrations (*Eqs. EA61
* and *
EA62
*).

#### Ca^{2+}* pump current.*

We included a sarcolemmal Ca^{2+} pump current (*I*
_{p,Ca}) formulation (*Eq.EA63
*) by following previous models (6, 22).

#### SR Ca^{2+}* handling.*

We implemented SR Ca^{2+} handling (*Eqs.EA64-EA73
*) using the two-compartment model (NSR and JSR) of Luo and Rudy (22) with a modification proposed by Friedman (12). The NSR subserves Ca^{2+} uptake, whereas the JSR governs release. As depicted in Fig. 1, the Ca^{2+}uptake current (*I*
_{up}) moves Ca^{2+} from the intracellular space to the network compartment, whereas the Ca^{2+} release current (*I*
_{rel}) releases Ca^{2+} from the junctional compartment to the myoplasm. A transfer current (*I*
_{tr}) moves Ca^{2+}from the NSR to the JSR. Ca^{2+} release from the junctional compartment is induced by Ca^{2+} flux into the myoplasm, with close coupling between sarcolemmal *I*
_{Ca} channels and *I*
_{rel} channels of the SR. The release current is thus a transient Ca^{2+} flux with a negative feedback that contributes to its rapid inactivation following a brief opening. Activation and inactivation gating for *I*
_{rel} were implemented as previously described (6).

#### Myoplasmic and SR Ca^{2+}* buffers.*

Ca^{2+} buffering is mediated by calmodulin and troponin in the myoplasm and by calsequestrin in the SR release compartment. Ca^{2+} binding by intracellular buffers was modeled as a dynamic process, as described in the Lindblad et al. (21) rabbit atrial myocyte model. Buffer concentrations and binding constants (see Table 1) were obtained from Luo and Rudy (22) and Rasmusson et al. (29).

## RESULTS

### Canine AP Model

The model-generated AP is shown in Fig.7. With the initial intracellular and extracellular ion concentrations given in Table2, repeated activation at 1 Hz results in stable intracellular concentrations of Na^{+} (11.8 mM), K^{+} (138.4 mM), Ca^{2+} (0.0001 mM), and Cl^{−} (29.3 mM) and a stable resting potential of −83.5 mV. Specific ionic currents are superimposed on the AP (Fig. 7, dotted line), and the [Ca^{2+}]_{i} transient is shown at the *lower right. *The fast, inward *I*
_{Na}is responsible for phase 0 upstroke. Inactivation of*I*
_{Na} and the onset of *I*
_{to},*I*
_{Kur,d}, and*I*
_{Cl,Ca} contribute to phase 1 repolarization, which lasts ∼10 ms, followed by a plateau established by a balance between inward *I*
_{Ca} and outward*I*
_{Kur,d}. The delayed rectifiers*I*
_{Kr} and *I*
_{Ks}, together with *I*
_{K1}, govern final repolarization.

### Variability of AP Morphology

Experimentally observed differences in AP morphology are often attributed to differences in the contributions of various ionic currents. We performed a sensitivity analysis to examine the relationship between AP morphology and densities of individual currents by varying the maximal conductance from 10 to 300% of the standard value in the model (Fig. 8). The control model AP is shown in each panel in bold. A 90% decrease in the conductance of a specified channel is represented by a dotted line, and a 300% increase in conductance is indicated by a dashed line (Fig. 8).

In general, AP duration (APD) increases when the conductance of a repolarizing current decreases (Fig. 8, *A–E*); however, the nature and degree of prolongation are different for each ionic current. For *I*
_{to}, the degree of APD prolongation saturates with decreasing conductance (Fig. 8
*A*), reaching a 17% prolongation of APD_{95} at a 25% conductance reduction. For *I*
_{Kur,d}, the plateau of the AP is raised as conductance decreases, accentuating spike and dome morphology and prolonging APD (Fig. 8
*B*). For both*I*
_{to} and *I*
_{Kur,d}, an increase in conductance progressively decreases APD and produces triangular APs. *I*
_{Kr} shows strong inward rectification positive to 0 mV, limiting its role in the early phases of the AP (Fig. 8
*C*); however, changes in*I*
_{Kr} produce important alterations in repolarization during late phase 2 and phase 3. Because of its slow activation and relatively small density, changes in*I*
_{Ks} do not alter APD_{95} by more than 5 ms in the model (data not shown). For *I*
_{K1}, small decreases in conductance (<50%) increase APD, with little change in the resting potential (Fig. 8
*D*). Larger decreases produce substantial depolarization, whereas increases in conductance reduce APD and cause very slight hyperpolarization. The contribution of*I*
_{Cl,Ca} to outward current during phase 1 is ∼5% of the combined densities of *I*
_{to} and*I*
_{Kur,d}; thus a 90% reduction in*I*
_{Cl,Ca} conductance causes only a 5-ms increase in APD_{95}, whereas a 300% increase in conductance decreases APD_{95} by ∼10 ms (Fig. 8
*E*). A reduction in*I*
_{Ca} progressively lowers the plateau and shortens APD, whereas an increase in conductance raises the plateau, accentuates spike and dome morphology, and prolongs APD (Fig.8
*F*).

### APD Adaptation to Rate

The adaptation of APD and refractoriness to changes in rate plays a significant role in the generation and maintenance of atrial arrhythmias (1, 10, 36, 38). Table3 lists model AP properties at rates between 0.1 and 5.0 Hz. As frequency increases, APD_{95}decreases from 338 to 130 ms. End-diastolic potential, AP amplitude, and overshoot are little affected. Figure9 shows model APs at 1, 2, and 3 Hz (Fig.9
*A*), which are in general agreement with experimental recordings (Fig. 9
*B*).

To identify the mechanism of rate adaptation in the model, we examined the effects of rate-dependent changes of individual selected currents on AP morphology. The variables characterizing a specific current at 4 Hz were substituted for those at 1 Hz, and a model AP was generated with the conditions for all other currents set as they were at 1 Hz. The resulting AP reflects the change that can be attributed to the rate-dependent response of the current under evaluation when the frequency is increased from 1 to 4 Hz, as illustrated in Fig.10. When the activation variables for*I*
_{Kr} and *I*
_{Ks} at 4 Hz are substituted into the model at 1 Hz, there is no significant AP change. Substitution of the voltage-dependent inactivation variable for*I*
_{Ca} (Fig. 10; denoted as *f* in Eq. 10) accounts for over 50% of the observed rate-dependent AP abbreviation. The addition of [Ca^{2+}]_{i}-dependent*I*
_{Ca} inactivation (by substitution of the [Ca^{2+}]_{i} values at 4 Hz) to voltage-dependent*I*
_{Ca} inactivation accounts virtually completely for AP rate adaptation.

### Regional Heterogeneity

AP heterogeneity plays an important role in the generation and maintenance of atrial reentrant arrhythmias (10, 26). Regionally determined differences in canine atrial AP morphology have been attributed to distinct variations in ionic current density (11); however, no critical analysis has been applied to determine whether experimentally observed regional differences in ionic currents are sufficient to explain the observed AP patterns. We therefore reproduced the experimentally observed differences in the model. Because model APs strongly resemble experimental recordings from pectinate muscle cells (11), we represented current densities from other regions in relationship to currents in pectinate muscle (Table 4) and performed simulations for each region with the use of ionic current density patterns recorded experimentally from that region. Figure11 shows the resulting model APs (Fig.11
*A*) along with digitally averaged experimental APs (Fig.11
*B*) obtained from all APs recorded from the corresponding regions by Feng et al. (11). Model AP morphologies are in good qualitative agreement with the experimental waveforms. Table5 compares APD_{95} in the model with the corresponding experimental averages in the four right atrial regions. The largest deviation of model APs from experimental recordings is in the appendage region, where the model AP had a less positive plateau and longer duration than the experimental average. This was the region with the greatest variation in experimentally recorded AP morphology (11), and overall morphology of the appendage AP in the model was much closer to the mean appendage experimental recording than to recordings from any other region.

### Tachycardia-Induced Atrial Ionic Remodeling

Rapid atrial activation provides a substrate for AF (13) that is characterized by reduced APD and APD rate adaptation (38). APD changes are accompanied by progressive decreases in *I*
_{Ca} and*I*
_{to} densities, summarized (relative to control values) in Table 6. To assess whether the reported changes in *I*
_{Ca} and*I*
_{to} are sufficient to explain the concomitant AP alterations, we reproduced the ionic alterations in the model.

Figure 12 shows the resulting model APs at 1 and 2 Hz (Fig. 12
*A*) and corresponding experimental recordings (*right*). The model shows APD abbreviation and loss of rate adaptation that is qualitatively similar to experimental observations. Quantitatively, rate-dependent APD adaptation was essentially abolished in myocytes of dogs subjected to rapid atrial pacing for 42 days, but still occurred in the model. The implementation of an additional reduction (to 10% of control values) of*I*
_{Ca} conductance (Fig. 12, *inset*) reproduced the experimentally observed loss of rate adaptation.

## DISCUSSION

We have developed a mathematical model of the canine atrial AP. The model produces APs that are consistent with experimental observations under a variety of conditions and permits critical analysis of the mechanisms and interactions of ionic currents that underlie AP properties in different situations.

### Behavior of the Model AP

#### Variability of AP morphology.

APs in isolated myocytes are known to exhibit a wide range of morphologies, even among cells isolated from small myocardial specimens. Wang et al. (35) identified three principal AP morphologies in atrial myocytes isolated from small human right atrial samples and attributed morphological variability to differences in the relative densities of repolarizing currents. Type 1 APs had a spike-and-dome waveform; type 2, an elevated, level plateau phase; and type 3, a triangular AP with little to no plateau region. Our sensitivity analysis (Fig. 8) demonstrates the effects of varying current densities on AP morphology and is consistent with the notion that the relative densities of currents may be responsible for differences in AP morphology.

#### Rate adaptation of the model AP.

Loss of atrial AP rate adaptation has been linked to electrical remodeling due to AF (13, 36) and vulnerability to AF induction in humans (1). Consistent with experimental observations (38), changes in *I*
_{Ca}were central to rate-dependent AP adaptation in the model (Fig. 10). Unlike our human atrial AP model (6),*I*
_{K} does not play a direct role in rate adaptation over 1–4 Hz in the present model. This discrepancy is likely due to the shorter and less positive plateau in the present model, because the amplitude of *I*
_{K} (particularly*I*
_{Kr}) is very sensitive to changes in the time the cell remains at voltages positive to the activation threshold (approximately −30 mV).

### AP Model Applications

#### Regional heterogeneity.

Heterogeneity of atrial refractoriness appears to play an important role in AF maintenance (10). Simulations of experimentally recorded regional differences in ionic conductances (11) produced APs with durations and morphologies closely resembling experimental observations (Fig. 11). Our model thus supports an ionic substrate as the basis for regional AP heterogeneity and constitutes the first quantitative effort of which we are aware to determine whether experimentally observed regional variations in ionic currents can account for regional AP differences.

#### AF-induced electrical remodeling.

The discovery of AF-induced remodeling (36) was an important advance in our understanding of AF mechanisms. Atrial tachycardia-induced ionic changes have been observed and suggested to account for the refractoriness alterations that are the hallmark of remodeling (38). Incorporating the ionic changes reported for remodeling in the model qualitatively reproduced experimental AP alterations (Fig. 12), indicating that the ionic alterations are likely an important contributor to AP remodeling. However, neither the overall AP abbreviation nor the reduction in AP rate adaptation was as striking in the model as those described experimentally. These findings suggest that the reported changes in *I*
_{Ca} and*I*
_{to} may be insufficient to account for tachycardia-induced changes in AP morphology or may be underestimated. Studies in dogs (14) and goats (31, 32) have suggested that electrical remodeling may be mediated by intracellular Ca^{2+} overload. Thus other ion transport processes, including the Na^{+}/Ca^{2+} exchanger (28), the Na^{+}-K^{+} pump, or the proton pump, which have not been studied experimentally in tachycardia-induced remodeling, may be involved. Inefficient coupling between sarcolemmal *I*
_{Ca} and SR Ca^{2+}release is associated with hypertrophy in rat ventricular myocytes (15, 30) and could be involved. Alternatively, the lack of full agreement with experimental findings may be due to inadequate expression by the model of changes in ionic processes at rapid rates.

### Comparison With Other Atrial AP Models

To our knowledge, this is the first canine-specific atrial AP model. In 1990, Rasmusson et al. (29) described a mathematical model of the bullfrog atrial cell. This model groups the voltage- and time-dependent repolarizing K^{+} currents into a single *I*
_{K}. The model (monophasic) AP displays a rounded peak and lacks a distinct phase 1. In addition, the model does not include a SR for intracellular uptake or release of Ca^{2+}.

Lindblad et al. (21) described a mathematical model of the rabbit atrial cell that incorporates intracellular Ca^{2+}handling by the SR and provides a more complete formulation of repolarizing K^{+} currents, including*I*
_{to}, *I*
_{Kr}, and*I*
_{Ks}. Unlike the present model,*I*
_{Ca} inactivation is unaffected by intracellular Ca^{2+}. The AP is very sensitive to rate changes over the range between 0.2 and 3 Hz because of the rate sensitivity of*I*
_{to} and its critical role in repolarization. Although that model faithfully reproduces properties of the rabbit atrial AP, species differences (particularly in the relative roles of*I*
_{to} and *I*
_{K}) limit its applicability to the dog.

Models of the human atrial AP have been published by Nygren et al. (27) and our group (6). Both models include similar ionic currents. Human atrial myocytes display a longer APD (by ∼60%) than in the dog. *I*
_{Ca} activates and inactivates faster in the dog model, reducing plateau voltage and duration, and consequently APD, in agreement with experimental AP recordings. The reduced plateau potential in the dog model also results in less contribution to repolarization from *I*
_{K}.

### Potential Limitations

Experimental ionic current estimates may be affected by bias that stems from the selection of cells with large currents and better viability, because cells with small currents may be considered inadequate for study and thus may not be represented in estimates of mean current densities. It is often necessary in formulating AP models to scale current densities to obtain appropriate AP properties (e.g., see Refs. 6, 17, and 27). In the present model,*I*
_{Ca} was the only current whose density required adjustment, with *I*
_{Ca} density in the model 30% of mean experimental values.

The mechanism of activation of *I*
_{Cl,Ca} is not yet fully understood. In the model, we formulated*I*
_{Cl,Ca} with a Ca^{2+}-dependent activation scheme that closely reproduced observed peak current-voltage relationships (Fig. 6). There is evidence that, in some experimental paradigms, reductions in *I*
_{Ca} are not accompanied by altered *I*
_{Cl,Ca} (38), suggesting that *I*
_{Cl,Ca} is regulated by additional factors not considered in the present model. Furthermore, the mechanisms governing intracellular Cl^{−} homeostasis are unclear. In the present model, we did not include a regulatory mechanism to counterbalance the outward movement of Cl^{−} via*I*
_{Cl,Ca}. Within the framework of the AP, however,*I*
_{Cl,Ca} is of limited importance (Figs. 7 and8
*E*), and the accumulation of intracellular Cl^{−}occurs very slowly at 1 Hz (∼1.6 mM/h or 6.43 mM in 4 h).

Alterations in *I*
_{Ks} produced little change in AP morphology or duration. Delayed rectifier currents are particularly sensitive to damage during cell isolation (39), so the limited role of *I*
_{Ks} in the model may be due to an underestimation of the current amplitude. Therefore, the finding that *I*
_{Kr} and *I*
_{Ks} play little role in AP rate adaptation should be interpreted with caution. Further experimental and theoretical work are needed to clarify better the properties and roles of *I*
_{Kr} and*I*
_{Ks} under physiological conditions.

### Potential Significance

The present model provides a useful tool for analyzing the role of ionic currents in canine atrial APs and arrhythmia mechanisms. Our model shows that interactions among currents play an important role in determining AP properties. Modulation of specific current conductances can have secondary effects on other currents via alterations in the voltage-time trajectory (Fig. 8). This is an important consideration when designing drugs that may alter ion channel properties.

AF is the most frequently encountered sustained arrhythmia in clinical practice. Key features associated with AF include loss of rate adaptation (1, 2) and increased heterogeneity of refractory periods (10, 26). Our model reproduces observed AP rate adaptation (Fig. 9) and thereby gives potential insights into underlying mechanisms. These investigations point to Ca^{2+}release from the SR, changes in intracellular [Ca^{2+}], and*I*
_{Ca} inactivation as important mediators of rate adaptation (Fig. 10) and support the role of an ionic substrate for regional AP heterogeneity as suggested by Feng et al. (11).

## Acknowledgments

We thank Lixia Yue, Danshi Li, and Jianlin Feng for providing experimental data and Drs. Normand Leblanc and Alain Vinet for helpful comments.

## MODEL FORMULATION

Some fractional equations require evaluation of a limit to determine their values at membrane potentials for which their denominator is zero Equation A1 Equation A2 Equation A3 Equation A4 Equation A5 Equation A6 Equation A7 Equation A8 Equation A9

### Equilibrium Potential

Equation A10

### Fast Na^{+}* Current*

Equation A11 Equation A12 Equation A13 Equation A14 Equation A15 Equation A17 Equation A18

### Time-Independent Inward Rectifier K^{+}* Current*

Equation A19

### Transient Outward K^{+}* Current*

Equation A20 Equation A21 Equation A22 Equation A23 Equation A24 Equation A25 Equation A26 Equation A27 Equation A28

### Ultrarapid Delayed Rectifier K^{+}*Current*

Equation A29 Equation A30 Equation A31 Equation A32 Equation A33 Equation A34 Equation A35 Equation A36 Equation A37 Equation A38

### Rapid Delayed Outward Rectifier K^{+}*Current*

Equation A39 Equation A40 Equation A41 Equation A42 Equation A43

### Slow Delayed Outward Rectifier K^{+}* Current*

Equation A44 Equation A45 Equation A46 Equation A47 Equation A48

### Sarcolemmal Ca^{2+}* Current*

Equation A49 Equation A50 Equation A51 Equation A52 Equation A53 Equation A54

### Ca^{2+}*-Activated Cl*^{−} Current

^{−}Current

Equation A55 Equation A56

### Na^{+}*-K*^{+}Pump Current

^{+}Pump Current

Equation A57 Equation A58 Equation A59

### Na^{+}*/Ca*^{2+}Exchanger Current

^{2+}Exchanger Current

### Background Currents

Equation A61 Equation A62

### Ca^{2+}* Pump Current*

Equation A63

### Ca^{2+}* Release Current From JSR*

Equation A64 Equation A65 Equation A66 Equation A67 Equation A68 Equation A69 Equation A70

### Transfer Current From NSR to JSR

Equation A71

### Ca^{2+}* Uptake Current by the NSR*

Equation A72

### Ca^{2+}* Leak Current by the NSR*

Equation A73

### Ca^{2+}* Buffers*

Equation A74 Equation A75 Equation A76

### Numerical Integration

At time step *p*, the updated value of a time-dependent variable is given by
Equation A77for χ = *V* or any time-dependent ionic concentration, and
Equation A78for *y* = any gating variable.

## Footnotes

This research was funded by grants from the Medical Research Council of Canada (S. Nattel) and the Natural Sciences and Engineering Research Council (M. Courtemanche). M. Courtemanche was supported by a Chercheur-Boursier Scholarship from the Fonds de Recherche en Santé du Québec. This work was included in the MScA thesis of R. J. Ramirez at the Université de Montréal.

Address for reprint requests and other correspondence: S. Nattel, Research Center, Montreal Heart Institute, 5000 Belanger St. E., Montreal, QC, Canada H1T 1C8 (E-mail:nattel{at}icm.umontreal.ca).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2000 the American Physiological Society