INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

THE ESTABLISHMENT of the original DiFrancesco-Noble (DN) model of cardiac myocyte electrophysiology spawned the development of numerous other dynamic models patterned after the DN formulation (5). These models explicitly include transmembrane ion channels and pumps, the intracellular calcium sequestering and release activity of the sarcoplasmic reticulum (SR), and changing intracellular ionic concentrations. The incorporation of concentration changes into the cardiac electrical model was an important development because such changes occur rapidly in the small volume of the cardiac cell and can have profound effects on action potential (AP) properties.

Dynamic models are finely tuned to reproduce electrophysiological behavior observed experimentally during observation periods of short duration. However, the original DN model was noted to produce a continual cycle-by-cycle (transient) accumulation (K⁺) or depletion (Na⁺, Ca²⁺) of intracellular ionic species (5). AP morphology and the pre- and postcycle values of the ionic concentrations were not constant as would be expected of steady-state behavior. Transient changes were negligible over a few cardiac cycles, but departures from the prestimulus initial concentrations were cumulative and became substantial over many cycles. Guan et al. (9) noted that the original DN equations are mathematically dependent (Jacobian is singular), and therefore that model fixed points are unstable. Varghese and Sell (32) showed that there exists a conservation principle latent in model equations that may permit certain mathematical features to cause computational difficulties and erroneous model performance. Although the abnormal behavior may occur only when parameters are outside the zone of physiological interest, it was argued that these features may introduce erroneous behavior and instabilities within the physiologically realizable range. Consequently, the validity of extended time simulations has been questioned (6, 26).

Time-dependent changes (transients) in ionic concentrations and APs following rate changes are well documented in the experimental literature (2, 7, 11-15, 27, 31, 38). Transients in mathematical models of the AP may represent artifacts of the system of equations used or may reflect physiological processes. We were unable to identify previous studies that characterized in detail time-dependent model transients and that compared model and experimentally observed AP transients. In the present study, electrophysiological transients in an ionic model of the canine atrial AP were studied with the following objectives: 1) to determine whether dynamic models reach stability during sustained pacing at a fixed rate; 2) to investigate the effects of stimulus current assignment; 3) to investigate the ionic basis of AP transients in the model; and 4) to compare experimental and model AP transients.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Model Implementation

The rate of change in the transmembrane potential (V) is given by

(1)

in time, where I_ion and I_stim are the total transmembrane ionic and stimulus currents, respectively, and C_m is the total membrane capacitance. Numerical integration was performed using a modified Euler method. After completion of the Cl transport formulation (described below), the total transmembrane ionic current is given by

(2)

I_ion includes contributions from the fast sodium current (I_Na), the inward rectifier (I_K1) and transient outward (I_to) potassium currents, the rapid (I_Kr) and slow (I_Ks) components of the classical delayed rectifier potassium current, L-type calcium current (I_Ca), a sarcolemmal calcium pump current (I_p,Ca), the Na⁺-K⁺-ATPase current (I_NaK), the Na⁺/Ca²⁺ exchanger (I_NaCa), and a calcium-activated chloride current (I_Cl,Ca), and background sodium (I_b,Na), calcium (I_b,Ca), and chloride (I_b,Cl) currents also contribute. The RNC AP most closely resembles right atrial pectinate muscle APs, because ionic current formulations were based on data from cells isolated from this region (16, 28). Shorter APs would be expected from a left atrial AP model, because I_Kr is typically larger in the left atria (17).

All simulations were performed with I in picoamps, V in millivolts, and C_m = 100 pF. A fixed time step of 5 and 20 µs was used in the presence and absence of stimulation, respectively. For reasons discussed in Stimulus Current Assignment (see RESULTS), all simulations were performed with the stimulus current attributed to potassium unless stated otherwise. All simulations were performed using double-precision arithmetic on Unix PC workstations.

Model Modification

The electroneutral Na⁺-Cl cotransporter was empirically formulated with a Hill function and is given by

(3)

where the conductance of the cotransporter (g_CT) = 0.115,  = E_Na  E_Cl (where E_Na is the equilibrium potential of Na⁺ and E_Cl is the equilibirium potential of Cl),  = 87.8251 and n = 4.

The background leakage current (I_b,Cl) is given by

(4)

where the conductance of the I_b,Cl (g_b,Cl) = 0.0018.

Because of the lack of kinetic data in the literature, g_b,Cl was approximated as twice the mean of the RNC background conductances. The cotransporter parameters were chosen to balance the magnitude of the leakage current at rest with rapid kinetics as intracellular pH is tightly regulated. The original RNC equations were based on short-term (several second) simulations. I_Ca was reduced to 30% of mean experimental values to obtain physiological AP durations (APDs). When we performed longer-term simulations, we found that the RNC parameters provided APDs that were too short, and that K⁺ currents had to be corrected [maximal I_to conductance (g_to) and I_Kur,d conductance (g_Kur,d) were reduced to 25 and 75% of original RNC values] to provide values closer to experimental observations. Conductance changes, CT_NaCl, and I_b,Cl were incorporated into the RNC equations to obtain the modified RNC (mRNC) model. Initial conditions for the mRNC model were obtained from the RNC initial conditions by allowing the revised equations to equilibrate for >5 min at rest (Table 1). Although the ionic mechanisms necessary for long-term stability were present following these modifications, it was noted that the equation singularity remained in the mRNC model. Therefore, the mRNC fixed points were also unstable and the equations still possessed the latent conservation principle as described.

Experimental Techniques for AP Recording

Preparations were stimulated with square-wave pulses (4 ms) delivered at 1.5 to 2 times diastolic threshold through bipolar Teflon-coated silver electrodes. The standard microelectrode techniques used in the present study have been described in detail elsewhere (24). Cellular membrane potentials were recorded from the endocardium using glass microelectrodes filled with 3 M KCl (8-20 M resistance) coupled to an Axoclamp 2B amplifier (Axon Instruments; Foster City, CA). Signals were converted into digital form by a Digidata 1200 series analog-to-digital converter (Axon Instruments) and displayed on a Pentium PC using Axotape version 2.0 software (Axon Instruments). Tissues were paced continuously at 400-, 300-, and 200-ms basic cycle lengths (BCLs) for 35 min at each BCL. AP recordings were obtained in each preparation from a minimum of five impalements 0-5, 15-20, and 30-35 min from the onset of stimulation at each BCL. A 5-min rest period followed stimulation at each BCL. During this time, preparations that did not beat spontaneously were paced at 1 Hz, and further control recordings were made to ensure the stability of the preparation. The 35-min pacing duration was chosen to allow measurement of the largest possible transients at three rates without compromising the viability of the tissues. Because results depended on the shape of the waveforms, precautions were taken to ensure that the stimulus current did not interfere. Each stimulus site was located at least several millimeters from the impaled cell (39). In addition, an interval of constant rest potential between the stimulus artifact and the recorded AP was confirmed (30). Light tension was applied to stabilize the preparations as needed. APD to 90% (APD₉₀) and 50% (APD₅₀) repolarization were measured with custom-made software and confirmed manually. Only recordings demonstrating <3% variation in interbeat end-diastolic potential were analyzed.

Statistical comparisons were performed on raw data with an exponential regression mixed-model analysis as described by Glantz and Slinker (8). Analysis of variance was applied for statistical comparison, with Bonferroni's correction used for post hoc tests. SAS release 6.12 (Cary, NC) was used to perform all statistical analyses. The level of statistical significance was set at P < 0.05.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Comparison of RNC and mRNC Ionic Transients

View larger version (23K): [in this window] [in a new window]   Fig. 1.   Steady-state model ionic transients. Transients in intracellular concentrations of K+ ([K+]i; A), Na+ ([Na+]i; B), Ca2+ ([Ca2+]i; C), and Cl ([Cl]i; D) are shown over 10 h of pacing at 2 Hz. Ramirez-Nattel-Courtemanche (RNC) concentrations failed to reach constant values. All modified RNC (mRNC) concentrations stabilized.

Although all mRNC concentrations appeared to be unchanging after 7 h of pacing, it was not clear whether the model system had reached absolute stability. To determine whether perfect stability was attained, numerical integration of all K⁺ currents and the AP waveform (measured to six significant figures) was performed after each hour of pacing for up to 10 h. All changes ceased between 6 and 7 h (not shown), confirming that model equations had reached absolute stability as suggested by the ionic concentrations (Fig. 1). In this way, the addition of Cl transporters was sufficient for the mRNC model to stabilize during sustained pacing, despite the inherent equation instabilities. The mRNC model was used for all subsequent simulations.

To ensure convergence of the integration scheme, the 10-h ionic transients were computed with two and four times reductions of the time step. No differences in [K⁺]_i were found between all time steps after 5 h. No differences between the 2.5- and 1.25-µs time steps were found for [Na⁺]_i after 5 h, and results with either differed from those with the 5-µs time step by <0.05%. [Ca²⁺]_i and [Cl]_i displayed slightly greater departures, but both converged toward the 5-µs solution. Overall, the 2.5- and 1.25-µs solutions differed from those with the 5-µs time step by <0.35% at all times, demonstrating that appreciable numerical error did not accumulate during the long pacing simulations and that the transients were indeed a property of the equation system.

Stimulus Current Assignment

View larger version (28K): [in this window] [in a new window]   Fig. 2.   Effects of stimulus current assignment on ionic transients. For each stimulus assignment (stm), resulting [K+]i (A), [Na+]i (B), [Ca2+]i (C), and [Cl]i (D) ionic transients are shown. All transients were altered by each assignment. Stimulus with K+ current assignment caused the smallest overall displacements.

Ionic Basis of [K+]_i Transient Reversal

View larger version (12K): [in this window] [in a new window]   Fig. 3.   Ionic mechanism of [K+]i transient reversal. After the onset of stimulation, [K+]i decreased (A) over many cycles because combined efflux currents (B) exceeded Na+-K+-ATPase current (INaK) influx. INaK gradually increased such that the direction of change in [K+]i reversed. Steady state was attained when influx and efflux were balanced. Action potential duration (APD) reduction (C) occurred very early and stabilized before [K+]i.

Unstable Fixed Points and Conservation Principle

To demonstrate the fixed-point instability, initial [K⁺]_i was decreased by 30, 60, and 90 mM, and the model was allowed to seek steady state by simulating 10 h without stimulation. If model fixed points were stable, concentrations would return to (or at least tend toward) the original unperturbed initial concentrations. The responses of all ionic species to the perturbations are shown in Fig. 4. It can be seen that the initial values of [K⁺]_i (Fig. 4A) were perturbed, whereas the profiles of [Na⁺]_i (Fig. 4B), [Ca²⁺]_i (Fig. 4C), and [Cl]_i (Fig. 4D) began at the same initial values. Each transient reached a new steady-state concentration (fixed point) with displacement proportional to the [K⁺]_i perturbation. The responses of [K⁺]_i and [Cl]_i were monotonic, whereas [Na⁺]_i and [Ca²⁺]_i changes reversed early in the simulation. This behavior resulted from the equation singularity and demonstrates the strict dependence of model solutions on initial conditions.

View larger version (25K): [in this window] [in a new window]   Fig. 4.   Unstable fixed points arising from mRNC equation singularity. Initial [K+]i was decreased by 30, 60, and 90 mM (A), and model equations were allowed to seek steady state by simulating 10 h without stimulation. At steady state, ionic transients were displaced in proportion to the [K+]i perturbations. 30K, 60K, and 90K indicate transients resulting from [K+]i decreases of 30, 60, and 90 mM, respectively.

The basic conservation principle (32) is given by

(5)

where u is the voltage, w is a concentration variable, and ₀ = C_m  1, ₁ = ₂ = (V_iF)  1, ₃ = (V_iF)  1, ₄ = (2V_iF)  1, ₅ = (2V_upF)  1, and ₆ = (2V_relF)  1.

The explicit form for the mRNC equations is given by

(6)

where V_i is the volume of the intracellular cytoplasm, and V_up and V_rel are the volumes of the SR uptake and release compartments, respectively. C₀ is a constant of integration.

This conservation principle explicitly states that the potential difference between the inside and outside of the cell is regulated by the flow of ions through the cell membrane. It is readily apparent from Eq. 6 that membrane potential may remain constant as long as the total intracellular charge is conserved, regardless of the ionic composition. In this way, the principle predicts that the unstable fixed points demonstrated by the displaced ionic transients in Fig. 4 would have appeared stable (i.e., transients would have returned to original unperturbed initial concentrations) if each [K⁺]_i decrease was offset with an equal but opposite increase in another cationic species.

This prediction was used to demonstrate the presence of the latent conservation principle in the mRNC equations. The initial value of [K⁺]_i was decreased by 30, 60, and 90 mM as before, with each charge deficit offset by increasing initial [Na⁺]_i by 30, 60, and 90 mM, respectively. The model was again allowed to seek steady state over 10 h without stimulation. The responses of all concentrations to the perturbations are shown in Fig. 5. It can be seen that [Ca²⁺]_i (Fig. 5C) and [Cl]_i (Fig. 5D) responded to the [K⁺]_i (Fig. 5A) and [Na⁺]_i (Fig. 5B) offsets with departures proportional to the displacement magnitudes. As expected, all transients returned completely to the original unperturbed initial conditions over 10 h (except for the response of [Cl]_i to the 90 mM perturbations, which returned to within >99% of the initial [Cl]_i value). The responses of [Na⁺]_i and [K⁺]_i were monotonic, whereas [Ca²⁺]_i and [Cl]_i reversed early in the simulation before returning toward baseline.

View larger version (26K): [in this window] [in a new window]   Fig. 5.   Conservation principle latent in the mRNC equations. Initial [K+]i was decreased by 30, 60, and 90 mM along with offsetting initial [Na+]i increases of 30, 60, and 90 mM, respectively. Model equations were allowed to seek steady state by simulating 10 h without stimulation as before. Although the perturbations caused large initial displacements of [K+]i (A), [Na+]i (B), [Ca2+]i (C), and [Cl]i (D), original initial concentrations were eventually recovered by all ionic species. 30K/Na, 60K/Na, and 90K/Na indicate transients resulting from initial [Na+]i and [K+]i displaced by 30, 60, and 90 mM, respectively.

Model AP Transients

View larger version (13K): [in this window] [in a new window]   Fig. 6.   Pacing-induced model AP transients with K+ stimulus assignment. APD at 90% (APD90; A) and and 50% repolarization (APD50; B) reductions are shown, with representative APs [300 ms basic cycle length (BCL)] from the 0-5 and 30-35 min pacing intervals (C). Marker at left indicates 0 mV.

View larger version (13K): [in this window] [in a new window]   Fig. 7.   Pacing-induced model AP transients with Na+ stimulus assignment. APD90 (A) and APD50 (B) reductions are shown, with representative APs (300 ms BCL) from the 0-5 and 30-35 min pacing intervals (C). Marker at left indicates 0 mV.

Ionic Basis of Model AP Transients

View larger version (28K): [in this window] [in a new window]   Fig. 8.   Ionic transients associated with time-dependent mRNC APD reductions during stimulation at cycle lengths indicated. Over the 35-min pacing interval, [K+]i (A) decreased, whereas [Na+]i (B), [Ca2+]i (C), and [Cl]i (D) increased.

The origin of the ionic transients was investigated. For the accumulation or depletion of ionic species the following inequality must hold

(7)

where i = K⁺, Na⁺, Ca²⁺, or Cl. I_i is the per-cycle sum of the charge carriers contributing to each concentration and is obtained by summing the relevant charge carriers from Eq. 2 with stoichiometric scaling of the pump and exchanger currents. For each ionic species, I_i is given by (with I_stim attributed to K⁺)

(8)

(9)

(10)

(11)

where SR(I_up, I_leak, I_rel) includes all Ca²⁺ fluxes between the SR and cytoplasm. Positive currents were defined as positive charge effluxes from the intracellular volume. The expressions for each charge carrier and the activity of the SR are fully described in Ref. 28.

To attain steady state at each stimulus rate, the cell attempts to remove the inequality in Eq. 7. Because certain charge carriers play only minor roles in the cardiac cycle, an equivalent contribution of each to I_i of Eq. 7 would not be expected. The currents that carried the greatest quantity of charge were most responsible for the magnitude of I_i and the resulting ionic transients. The greatest changes in current magnitude toward this end contributed most to neutralizing the concentration drift. To determine the total charge carried by all membrane currents and the functional changes in current magnitude over the pacing interval, numerical time integration of each current in Eqs. 8-11 was performed over one 300-ms cycle after 2.5 and 32.5 min of pacing. Integration results and the corresponding I_i are given in Tables 2-5. In the tables,  indicates the magnitude of the contribution of functional changes to ionic stabilization. I_i indicates the overall charge imbalance at each time point, and I_i indicates the magnitude of combined functional changes toward ionic stabilization.

K+ and 35-min transients. All K⁺ currents (Table 2) adjusted toward a new [K⁺]_i steady state, except the stimulus current, which has a fixed magnitude in the RNC model. The largest K⁺ current was efflux by I_K1 after 2.5 and 32.5 min. Over the pacing interval, reduced efflux by I_Kur,d contributed most toward [K⁺]_i stabilization, although reduced effluxes by I_K1 and I_to were also important to limiting [K⁺]_i depletion. In Table 2, negative I for membrane currents indicates that current magnitude decreased during the pacing interval. Positive I_K indicates increased conservation of intracellular K⁺.

Na+ and 35-min transients. All Na⁺ charge carriers (Table 3) adjusted toward a new [Na⁺]_i steady state, with the exception of I_b,Na. The largest Na⁺ charge carrier at both times was efflux by I_NaK. The greatest contributor to the pacing-induced increase in [Na⁺]_i was I_NaCa, and functional changes in the balance of reverse-mode I_NaCa exchange also contributed most toward [Na⁺]_i stabilization. In Table 3, negative I_Na indicates decreased accumulation of intracellular Na⁺.

Ca²+ and 35-min transients. For Ca²⁺ (Table 4), I_Ca and I_p,Ca adjusted toward a new [Ca²⁺]_i steady state, whereas I_NaCa and I_b,Ca opposed it. The overall magnitude of I_NaCa preferentially attenuated in favor of establishing a new [Na⁺]_i steady state. The largest Ca²⁺ current was efflux by I_NaCa after 2.5 and 32.5 min. Reduced influx by I_Ca contributed most toward establishing a new [Ca²⁺]_i steady state. The immediate drop following the onset of stimulation was due to large Ca²⁺ efflux carried by I_NaCa. The preceding train of conditioning pulses at 1 Hz before the onset of rapid stimulation reduced this early redistribution. The greater final magnitude of [Ca²⁺]_i at shorter BCLs evinced the rate dependence of Ca²⁺ loading. SR fluxes were balanced during each cycle, and a net SR flux did not contribute appreciably to the Ca²⁺ overload. In Table 4, negative I_Ca indicates decreased accumulation of intracellular Ca²⁺.

Cl and 35-min transients. Finally, the magnitudes of all Cl charge carriers (Table 5) adjusted toward establishing a new [Cl]_i steady state, except I_Cl,Ca, which was constrained to increase with increasing [Ca²⁺]_i, thereby contributing to the delayed [Cl]_i steady state seen in Fig. 8. The largest Cl charge carriers were influx by CT_NaCl at both times. Increased efflux by I_b,Cl also contributed most toward establishing a new [Cl]_i steady state, although decreased influx by CT_NaCl was only slightly less important. Compared with other functional changes contributing to charge balance, adjustment of Cl charge carriers were of minor importance. In Table 5, negative I_Cl indicates decreased accumulation of intracellular Cl.

Contribution of concentration transients to AP reduction. To determine the relative contribution of each 35-min concentration transient (Fig. 8) to pacing-induced AP reduction, the 2.5-min AP was computed after clamping initial concentrations with values sampled 2 ms before the AP upstroke after 32.5 min of pacing (BCL 300 ms). Results are shown in Fig. 9. The addition of rapid pacing-induced changes in [K⁺]_i (1.5% decrease) prolonged the 2.5-min measurement of APD₉₀ by 1 ms (1.0% prolongation) but did not change APD₅₀ (Fig. 9A). Pacing-induced changes in [Na⁺]_i (18% increase) accounted for 100% and 121% of the overall APD₅₀ and APD₉₀ reductions, respectively (Fig. 9B). The reduction of APD₉₀ was beyond that of the 35-min AP because these effects were augmented in the absence of the APD-preserving [K⁺]_i decrease and [Ca²⁺]_i change. Because [Ca²⁺]_i is not constant over the AP and changes in Ca²⁺ handling strongly affect the [Ca²⁺]_i time course, myoplasmic [Ca²⁺]_i changes (24.7% increase) and changes in SR uptake (13.5% increase) and release (18.1% increase) compartments were also included. APD₅₀ was unchanged, whereas APD₉₀ was prolonged by 3 ms (3.1% prolongation). Pacing-induced changes in [Cl]_i (13.5% increase) did not alter APD. It was thus concluded that pacing-induced increased [Na⁺]_i contributes significantly to AP shortening, whereas changes in [K⁺]_i and Ca²⁺ handling oppose this effect. This result is consistent with Fig. 7, where the Na⁺ stimulus assignment markedly shortened APD relative to the K⁺ assignment, because of greater cell intracellular Na⁺ loading with Na⁺ as the stimulus current and a consequent outward shift in I_NaCa.

View larger version (12K): [in this window] [in a new window]   Fig. 9.   Effects of rate-related ionic concentration changes on APs. Results were obtained by substituting 32.5-min values for intracellular concentration of each species (K+ in A, Na+ in B, Ca2+ in C, and Cl in D) into the equations for the AP at 2.5 min. For Ca2+, 32.5-min Ca2+-handling parameters were also substituted, because intracellular Ca2+ during the AP is determined largely by Ca2+ handling in the sarcoplasmic reticulum.

Contribution of functional current changes to AP reduction. To determine the relative contribution of the pacing-induced I_Ca,L decrease to AP shortening, it was first necessary to determine the changes in calcium handling at 32.5 min without the influence of the abbreviated 32.5-min AP waveform. This was accomplished by computing the 32.5-min calcium transient ([Ca²⁺]_i changes during single AP) with clamping of the 2.5-min AP waveform. I_Ca,L underlying the 2.5-min AP with clamping of the modified calcium transient (for I_Ca,L only) is shown in Fig. 10A. The peak magnitude of the current was slightly reduced due to increased [Ca²⁺]_i-induced inactivation of I_Ca,L, a mechanism consistent with experimental studies (20, 21). These changes were not sufficient to alter APD₅₀ or APD₉₀ of the 2.5-min AP (Fig. 11A). To determine the relative contribution of I_NaCa changes to the pacing-induced AP reduction, the 2.5-min exchanger current was computed with clamping of the modified calcium transient (used for I_Ca,L above) and 32.5-min [Na⁺]_i, both for I_NaCa only, producing the current and AP changes illustrated in Figs. 10B and 11B, respectively. The 2.5-min AP was significantly shortened by the I_NaCa changes occurring at 32.5 min, accounting for 83% and 104% of APD₉₀ and APD₅₀ reduction, respectively. Clearly, I_NaCa changes contributed importantly to rate-dependent APD alterations. To determine the relative contribution of the pacing-induced I_NaK increase to the corresponding AP reduction (BCL 300 ms), the 2.5-min AP was computed with the 32.5-min values of [K⁺]_i and [Na⁺]_i clamped for I_NaK only. The altered I_NaK waveform and resulting AP changes are shown in Figs. 10C and 11C, respectively. APD₅₀ of the 2.5-min AP was unchanged, whereas the increased net outward current (also Tables 2 and 3) accounted for 4.2% of APD₉₀ reduction. Analysis of other currents showed that they changed minimally between 2.5 and 32.5 min and did not contribute significantly to AP alterations over this interval.

View larger version (9K): [in this window] [in a new window]   Fig. 10.   Effects of current parameters at 32.5 min on current during 2.5 min AP. A: ICa,L; B: INaCa; C: INaK. For each current, the defining parameters at 32.5 min were substituted into the 2.5-min AP equations, indicating the effects of pacing-induced parameter changes on current during the AP.

View larger version (8K): [in this window] [in a new window]   Fig. 11.   Effects of pacing-induced changes in each current on the AP. Currents obtained in Fig. 10 were substituted individually into the AP at 2.5 min. A: ICa,L clamp; B: INaCa clamp; C: INaK clamp. Resulting AP indicates the effect that pacing at 32.5 min had on the AP via the current indicated. Results show that INaCa played the largest role, in itself accounting for most of the rate-related AP changes.

Experimental and Model AP Transients

View larger version (16K): [in this window] [in a new window]   Fig. 12.   Pacing-induced experimental AP transients. APD90 (A) and APD50 (B) reductions are shown with representative APs (300 ms BCL) from the midpoint of 0-5 and 30-35 min pacing intervals (C). Marker at left indicates 20 mV.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

We have shown that dynamic models with unstable fixed points may stabilize, provided complete homeostatic mechanisms for all ionic species are present, including the stimulus charge. Model stability depended on the assignment of the stimulus current, and K⁺ stimulus assignment produced the smallest perturbations in the equilibrium of all ionic species. Model AP transients agreed qualitatively with experiments, and therefore were not appreciably influenced by potential model artifacts due to unstable fixed points. AP reduction was due to changes in the Na⁺/Ca²⁺ exchanger function related to accumulation of intracellular Na⁺. APD reduction was opposed by transient changes in [K⁺]_i and [Ca²⁺]_i.

Comparison With Previous Cardiac Models

Ionic drift is a common property of ionic models based on the original DN formulation scheme (5). Similar models include the human atrial myocyte models by Nygren et al. (26) and Courtemanche et al. (3), rabbit atrial myocyte model by Lindblad et al. (18), rabbit sinoatrial node model by Demir et al. (4), the ventricular myocyte models of Nordin (25) and the phase II Luo and Rudy (LR2) model (19), and the bullfrog atrial myocyte model by Rasmusson et al. (29). Models by Jafri et al. (10) and Winslow et al. (37) are also similar because they represent the Luo and Rudy equations with modified subcellular Ca²⁺ handling.

None of these earlier models monitored [Cl]_i. Lindblad et al. (18) included a formulation for I_b,Cl, but [Cl]_i was held constant such that Cl fluxes did not disturb the ionic equilibrium. Nygren et al. (26) added a small electroneutral flux of Na⁺ to achieve long-term ionic homeostasis. It was suggested that this flux could be accounted for by a cotransport system with Cl similar to the one implemented in the present study; however, these mechanisms were not modeled. Because these other models account for only [K⁺]_i, [Na⁺]_i, and [Ca²⁺]_i, ionic equilibrium may theoretically be attained by I_NaK and I_NaCa. Rather than fixing [Cl]_i, we chose to provide for [Cl]_i equilibrium by formulating an empirical model of Cl transport based on physiological mechanisms. This allowed changes in I_Cl,Ca to influence the ionic transients.

Transients in these earlier models were qualitatively similar to those of the mRNC model. In the model by Lindblad et al. (18), a small increase in [Na⁺]_i was noted over each pacing cycle, and a slow downward drift of [Na⁺]_i followed the termination of pacing. Rasmusson et al. (29) reported that [Na⁺]_i decreased over 2 min in the absence of stimulation. [Na⁺]_i increased and [K⁺]_i decreased over 40 s of pacing, and the magnitude of the transients was greater at faster rates. In the LR2 model, [K⁺]_i was reported to decrease, whereas [Ca²⁺]_i and [Na⁺]_i increased during 60 s of pacing at 2 Hz (19). The transients were used to illustrate the redistribution of ions that may result from overdrive and cause suppression of pacemaker cells. In the model by Courtemanche et al. (3), it was noted that the ionic balance was disturbed by periodic stimulation and that slow changes in intracellular ionic concentrations still persisted after several minutes. Initial transients were attributed to kinetic rate adaptation of the currents and dissipated over a few seconds. The small subsequent slow changes in ionic concentrations resulted in slow changes in the shape of the AP. To ensure that the AP morphology took into account kinetic adaptation at the specified frequency but did not include the effects of concentration changes, simulation results were presented after 12 s of pacing from rest. Similarly, Demir et al. (4) showed results after simulation of at least eight cycles. To avoid the effects of concentration changes on model results, Nygren et al. (26) added a small electroneutral inward flux of Na⁺ (discussed above), and Dokos et al. (6) modified the maximum I_NaK and time constant for uptake of intracellular calcium of the original DN equations.

In the model by Nordin (25), [Na⁺]_i was monitored over 50 min of sustained pacing. Pacing was initiated at 1 Hz for 10 min and was then increased without interruption to 2 and 4 Hz for 10 min at each rate. [Na⁺]_i markedly increased and then began to plateau during each interval. At 30 min, pacing was slowed to 2 Hz and then 1 Hz for 10 min at each rate. [Na⁺]_i decreased along an exponential time course as the stimulus rate was slowed and returned to baseline after stimulation was stopped. The Nordin model also generated a fourfold increase in the size of myoplasmic [Ca²⁺]_i transients as the stimulation rate increased and diastolic [Ca²⁺]_i more than doubled. The role of the ionic transients on AP shortening was also studied. As in the mRNC model, [Na⁺]_i-independent AP shortening over several minutes of stimulation at rapid rates was caused in part by a reduction of Ca²⁺ current secondary to increased myoplasmic [Ca²⁺]_i (25). AP shortening was largely caused by an outward shift of I_NaCa as [Na⁺]_i increased (25).

Comparison With Experimental Findings

Stimulus Current Assignment

There is no experimental evidence on which to allocate the stimulus current carrier for model simulations. The most likely charge carriers are K⁺ and Na⁺, which are the cations that are by far the most concentrated in the intracellular and extracellular compartments, respectively. The resting ionic concentrations in the heart favor Na⁺ as an inward current carrier. However, during impulse propagation, cell-to-cell coupling is effected by pore-forming gap junctions, which are highly permeable to intracellular cations. Because K⁺ is by far the most concentrated mobile intracellular ionic species, it is likely that K⁺ movement plays an important role in excitatory cell-to-cell transmission. In addition, when a portion of the membrane is depolarized, the adjacent membrane may respond by depolarization via local movement of subsarcolemmal K⁺. The results with Na⁺ as the charge carrier (Fig. 7) resulted in unphysiologically short action potentials, because Na⁺ loading during activity favored reverse-mode outward I_NaCa (Figs. 9-11). Because attributing the stimulus current to K⁺ reproduced experimental observations of [K⁺]_i changes (Fig. 2), because the K⁺ assignment best preserved the intracellular ionic equilibrium (Fig. 2), and because AP transients were most physiological for the K⁺ stimulus carrier (Fig. 6), K⁺ stimulus attribution was the preferred assignment. Assignment of stimulus to an ionic species is necessary for stability. Physiologically, the stimulus current may actually be carried by more than one ionic species; however, in the absence of knowledge about the relative participation of different ionic species, it is reasonable and practical to attribute stimulus current to a single ionic species for long-term simulations. It is possible that lack of stability in related models may be corrected by the appropriate stimulus current assignment.

Novel Aspects and Potential Significance

Sustained rapid pacing is known to induce electrophysiological remodeling, mimicking chronic atrial fibrillation (AF) (22, 35). Tachycardia-induced remodeling shortens APD and abolishes AP rate adaptation, thereby enhancing atrial arrhythmogenicity ("AF begets AF") (35, 36). In atrial myocytes of dogs atrially paced at 400 beats/min for 6 wk, Yue et al. (40) found that the densities of I_Ca and I_to were progressively decreased by 69 and 65%, respectively, as a result of downregulation of mRNA encoding I_Ca and I_to -subunits (41). The most extensive remodeling occurred early and gradually attenuated toward the end of the pacing period. In contrast, no remodeling of I_Kr, I_Ks, I_K1, I_Kur,d, or I_Cl,Ca was found. The present study provides insights into the functional events following the onset of rapid pacing. Results here indicate that ionic transients cause short-term APD reductions (over 30 min) that precede remodeling per se. These functional changes may contribute to early AF-promoting effects of atrial tachycardia (23).

Potential Limitations

Agreement between model and experiment was theoretically limited. The explicit form of the conservation principle (Eq. 5) includes a constant of integration, C₀. The physical significance of C₀ may be understood by recognizing that the conservation principle is Faraday's principle
where Q₀ is extracellular charge and Q_i is intracellular charge rewritten in terms of ionic concentration. Here, Q_o is fixed, and Q_i changes throughout the cardiac cycle. Although the model accounts for the concentrations of the ionic species predominantly responsible for the AP, other ionic species (e.g., HCO, SO, PO, and Mg²⁺) are found within the cardiac myocyte as well as extracellularly. C₀ contains the combined effect of additional ionic species and may be expressed as

(12)

where Q₀ and Q_i indicate extracellular and intracellular charge, respectively. The presence of these additional ionic species in physiological systems necessarily implies discrepancies with model results.

It was recognized, as discussed above, that experimental AP recordings were made from a tissue syncytium, whereas the model represents an isolated cell. The role of one-dimensional propagation effects was investigated by modifying Eq. 1 to obtain the reaction-diffusion system

(13)

in time and one space dimension, where D = 0.001 cm²/ms and ² is the second-derivative Laplacian operator in space.

The finite difference equation resulting from Eq. 13 is

(14)

where V is the voltage at time step n and location j. t and x are the time and space discretization steps fixed at 5 µs and 0.025 cm, respectively. I_couple is the coupling current between adjacent cells with a 3-point central difference approximation to the Laplacian given by

(15)

Preliminary cable simulations with I_couple attributed to K⁺ revealed that transients were augmented in coupled segments, and that the net charge carried by I_couple was nearly 30 times less than I_stim. The finding that reduced net stimulus current augmented the ionic transients is consistent with the exaggerated transients observed in the single cell when the stimulus current was not assigned (Fig. 2). Although it may have been meaningful to study model stability and fully characterize model transients in the context of propagation, this was not done for reasons of computational tractability. Wall time required to pace a 50-cell cable for 1 h was on the order of weeks, and simulations of this type may not be parallelized. It would be expected that transients in a propagated model would also reach absolute stability, although substantially more time would likely be required in a distributed system.

We assumed that extracellular ionic concentrations were constant. This is a reasonable assumption for isolated cells. However, the type of experimental preparation may influence ionic transients (13), and in a multicellular preparation ion accumulation/depletion phenomena can occur and contribute to rate-related changes. Also, biochemical and molecular changes may occur in vivo that are not reproduced in the purely ionic mathematical model. Finally, the experimental recording method also includes potential artifacts, including possible tissue ischemia during perfusion with crystalloid solutions that have limited oxygen-carrying capacity and tissue damage during isolation and experimental preparation.