## Abstract

Atrial fibrillation, a common cardiac arrhythmia, is promoted by atrial dilatation. Acute atrial dilatation may play a role in atrial arrhythmogenesis through mechanoelectric feedback. In experimental studies, conduction slowing and block have been observed in acutely dilated atria. In the present study, the influence of the stretch-activated current (*I*_{sac}) on impulse propagation is investigated by means of computer simulations. Homogeneous and inhomogeneous atrial tissues are modeled by cardiac fibers composed of segments that are electrically and mechanically coupled. Active force is related to free Ca^{2+} concentration and sarcomere length. Simulations of homogeneous and inhomogeneous cardiac fibers have been performed to quantify the relation between conduction velocity and *I*_{sac} under stretch. In our model, conduction slowing and block are related to the amount of stretch and are enhanced by contraction of early-activated segments. Conduction block can be unidirectional in an inhomogeneous fiber and is promoted by a shorter stimulation interval. Slowing of conduction is explained by inactivation of Na^{+} channels and a lower maximum upstroke velocity due to a depolarized resting membrane potential. Conduction block at shorter stimulation intervals is explained by a longer effective refractory period under stretch. Our observations are in agreement with experimental results and explain the large differences in intra-atrial conduction, as well as the increased inducibility of atrial fibrillation in acutely dilated atria.

- atrial fibrillation
- excitation-contraction coupling
- impulse propagation
- stretch-activated current

atrial fibrillation (AF) is a common cardiac arrhythmia (33). An important risk factor for AF is chronic atrial dilatation (38, 56), whereas experimental studies indicate a role of acute atrial dilatation in the initiation of atrial arrhythmia (2, 34, 39, 48, 50). Conduction slowing and shortening of the refractory period in acutely dilated atria have been reported (6, 16, 40). Eijsbouts et al. (10, 11) found, in addition to conduction slowing, an increased occurrence of intra-atrial block. Hu and Sachs (15) and Kohl and Sachs (28) hypothesize that stretch-induced changes in electrophysiological behavior can be explained by stretch-activated channels (SACs). In the present simulation study, we investigate this hypothesis for atrial impulse propagation.

Several models have been proposed to describe SACs based on experimental observations (12, 45, 52, 58, 59). Similar models have been applied in large-scale computer simulations to investigate the effect of stretch on defibrillation (53) and the termination of ventricular tachycardia by means of precordial thump (31). Since cardiomechanics are not considered in these studies, the stretch-activated current (*I*_{sac}) is not influenced by contraction. Models of the ventricles in which contraction is triggered by electrical activation describe stimulation from the Purkinje system (24, 55) and epicardial stimulation (25, 54). In these studies, mechanical deformation is triggered by electrical activation. However, mechanoelectric feedback, i.e., the effect of mechanical deformation on the electrophysiology, is not considered. To investigate the influence of mechanical deformation on impulse propagation, a strong coupling between cardiomechanics and electrophysiology is required, as proposed elsewhere (32, 35, 36). In these studies, tissue conductivity is directly affected by mechanical deformation, and the amount of *I*_{sac} is related to local deformation of the cardiac tissue. Physiological details, such as ionic membrane currents, intracellular Ca^{2+} handling, and cross-bridge formation, are not considered in these models.

In the present study, we investigate the role of *I*_{sac} in conduction slowing and block as observed in acutely dilated atria. We apply a discrete bidomain model with strong coupling between cardiomechanics and cardiac electrophysiology. Our model describes ionic membrane currents, Ca^{2+} storage and release from the sarcoplasmic reticulum (SR), and cross-bridge formation. In contrast to all other multicellular models, contractile forces are directly coupled to free Ca^{2+} concentration, as well as sarcomere length. In our model, the amount of *I*_{sac} is related to local stretch and may change during contraction. We performed simulations of homogeneous and inhomogeneous cardiac fibers under stretch to quantify the conduction velocity in the presence of *I*_{sac}. We observed conduction slowing, a longer effective refractory period (ERP), and (unidirectional) conduction block with increasing stretch. Furthermore, we found that contraction of early-activated fiber segments can lead to conduction block in later-activated segments. The observed phenomena are in agreement with experimental observations and provide an explanation for the increased inducibility of AF in acutely dilated atria.

## METHODS

In the present study, we apply our discrete bidomain model, the cellular bidomain model (29, 30), which describes active membrane behavior, as well as intercellular coupling and interstitial currents, and has been extended to model cardiac tissue mechanics and *I*_{sac}. We describe the extensions to our model of cardiac electrophysiology, in particular the influence of stretch on fiber conductivity, our model of the *I*_{sac}, the Ca^{2+}-force relation, the mechanical behavior of a single segment, and the mechanical behavior of a cardiac fiber. Furthermore, the numerical integration scheme is described, and an overview of the simulations is given.

### Modeling Cardiac Electrophysiology

In the cellular bidomain model, the cardiac tissue is subdivided into segments, each with its own membrane model describing the ionic membrane currents (29, 30). The state of each segment is defined by the intracellular potential (*V*_{int}), the extracellular potential (*V*_{ext}), and the state of the cell membrane, which is expressed in gating variables and ion concentrations. The membrane potential (*V*_{mem}) is defined by (1) Intracellular and extracellular currents between adjacent segments are related to intracellular and extracellular conductivities and satisfy Ohm's law. Exchange of current between the intracellular and extracellular space occurs as transmembrane current (*I*_{trans}) (2) where χ is the surface-to-volume ratio needed to convert *I*_{trans} per unit membrane surface to *I*_{trans} per unit tissue volume, *C*_{mem} is the membrane capacitance, which is typically 1 μF/cm^{2} for biological membranes (1), and *I*_{ion} is ionic current (expressed in μA/cm^{2} membrane surface or pA/pF when *C*_{mem} = 1 μF/cm^{2}). *I*_{ion} depends on *V*_{mem}, gating variables, and ion concentrations (see below). Impulse propagation is related to the longitudinal conductivity parameters *g*_{int} and *g*_{ext}. The bidomain parameters used for the present study are from Henriquez (13) and are based on measurements by Clerc (7) (Table 1).

To incorporate *I*_{sac}, we modified the model of the human atrial action potential (AP) of Courtemanche et al. (8). The total ionic current is given by (3) where *I*_{Na} is fast inward Na^{+} current, *I*_{K1} is inward rectifier K^{+} current, *I*_{to} is transient outward K^{+} current, *I*_{Kur} is ultrarapid delayed rectifier K^{+} current, *I*_{Kr} is rapid delayed rectifier K^{+} current, *I*_{Ks} is slow delayed rectifier K^{+} current, *I*_{Ca,L} is L-type Ca^{2+} current, *I*_{p,Ca} is Ca^{2+} pump current, *I*_{NaK} is Na^{+}-K^{+} pump current, *I*_{NaCa} is Na^{+}/Ca^{2+} exchanger current, and *I*_{b,Na} and *I*_{b,Ca} are background Na^{+} and Ca^{2+} currents (8). The model for the ionic and pump currents, including handling of intracellular Ca^{2+} concentration ([Ca^{2+}]_{i}) by the SR, is adopted from the model of Courtemanche et al. (8).

### Influence of Stretch on Fiber Conductivity

The intracellular and extracellular conductivities (*g*_{int} and *g*_{ext}) may change during stretch or contraction of the fiber. Under stretch, the length of the cells increases and the cross-sectional area decreases, leading to a reduced fiber conductivity. To quantify the changes in *g*_{int} and *g*_{ext} (mS/cm), we assume that the resistivity of the intracellular space (*R*_{int} = 1/*g*_{int}, Ω·cm) is determined partly by the myoplasmic resistivity (*R*_{myo}) and partly by the gap-junctional resistivity (*R*_{junc}) (4) For the nonstretched fiber, we define *g*_{int} = *g*_{int 0} = 1/*R*_{int 0}, *R*_{myo} = *R*_{myo 0}, and *R*_{junc} = *R*_{junc 0}. When the fiber is stretched with stretch ratio λ, cell length increases and cross-sectional area decreases (assuming that cell volume is conserved). Since *R*_{myo} is proportional to the length and inversely proportional to the cross-sectional area of the cell, we obtain (5) On the basis of the assumption that the total number of gap junctions in the fiber does not change under stretch, the number of gap junctions per length unit decreases proportionally with λ, which leads to (6) If *Eqs. 5* and *6* are combined, *g*_{int} is related to λ by (7) For the extracellular domain, we assume that *g*_{ext} is related to *g*_{ext 0} and λ by (8) Chapman and Fry (5) determined that 52% of the total resistivity was attributed to gap-junctional resistance in frog myocardial cells (*R*_{junc 0}/*R*_{int 0} = 0.52). Since these cells are longer (131 μm) (5) than human atrial cells (94 μm) (37), we estimate that *R*_{junc 0}/*R*_{int 0} = 0.6 for human atrial myocardium (Table 1).

### I_{sac}

On the basis of experimental observations, we assume that *I*_{sac} in atrial myocytes is a nonselective cation current with a near-linear current-voltage relation (26). The reversal potential is −3.2 mV for rat atrial myocytes (26). In our model, *I*_{sac} is permeable to Na^{+}, K^{+}, and Ca^{2+} and is defined by (9) where *I*_{sac,Na}, *I*_{sac,K}, and *I*_{sac,Ca} represent the Na^{+}, K^{+}, and Ca^{2+} contributions, respectively, to *I*_{sac}. These currents are defined by the constant-field Goldman-Hodgkin-Katz current equation (22) (10) (11) (12) where *P*_{Na}, *P*_{K}, and *P*_{Ca} denote the relative permeabilities to Na^{+}, K^{+}, and Ca^{2+}, *z*_{Na}, *z*_{K}, and *z*_{Ca} represent the ion valences, and *F* is Faraday's constant, *R* is the universal gas constant, and *T* is temperature (310 K) (8).

The conductance (*g*_{sac}) depends on λ as follows (13) where *G*_{sac} is the maximum membrane conductance, *K*_{sac} is a parameter to define the amount of current when the cell is not stretched [λ = 1, sarcomere length (*l*_{s}) = 1.78 μm], and α_{sac} is a parameter to describe the sensitivity to stretch. *K*_{sac} and α_{sac} are from Zabel et al. (58) (Table 1).

The reversal potential (*E*_{sac}) can be obtained by solving the following equation for *V*_{mem}: *I*_{sac,Na} + *I*_{sac,K} + *I*_{sac,Ca} = 0. In the present study, we consider two cases: *P*_{Na}:*P*_{K}:*P*_{Ca} = 1:1:1, with *E*_{sac} = −0.2 mV, and *P*_{Na}:*P*_{K}:*P*_{Ca} = 1:1:0, with *E*_{sac} = −0.9 mV. In both cases, *I*_{sac} has a near-linear current-voltage relation (Fig. 1).

To describe the influence of *I*_{sac,Na}, *I*_{sac,K}, and *I*_{sac,Ca} on intracellular Na^{+}, K^{+}, and Ca^{2+} concentrations ([Na^{+}]_{i}, [K^{+}]_{i} and [Ca^{2+}]_{i}), respectively, we replace *Eqs. 21*–*25* of the model of Courtemanche et al. (8) by (14) (15) (16) (17) (18) where *C*_{m} is the membrane capacitance of a single atrial myocyte (100 pF) (8), *F* is Faraday's constant, V_{i} is the intracellular volume (13,668 μm^{3}) (8), V_{up} and V_{rel} are the volumes of the SR uptake and release compartments, respectively, *I*_{up,leak}, *I*_{up}, and *I*_{rel} represent the SR currents, [Trpn] is troponin concentration, [Cmdn] is calmodulin concentration, and *K*_{m} is the half-saturation constant. *Equation 18* is *Eq. 25* in the model of Courtemanche et al. (8) and represents the influence of Ca^{2+} buffering in the cytoplasm mediated by troponin ([Ca^{2+}]_{Trpn}) and calmodulin ([Ca^{2+}]_{Cmdn}) on [Ca^{2+}]_{i}.

### Modeling the Ca^{2+}-Force Relation

Rice et al. (42–44) proposed five models of isometric force generation in cardiac myofilaments. To model the Ca^{2+}-force relation in the present study, we apply their *model 4*, which is based on a functional unit of troponin, tropomyosin, and actin. The binding of Ca^{2+} to troponin is described by two states: unbound troponin and Ca^{2+} bound to troponin. Tropomyosin can be in one of six states: nonpermissive with 0 and 1 cross bridges (*N0* and *N1*) and permissive with 0, 1, 2, and 3 cross bridges (*P0*, *P1*, *P2*, and *P3*). The permissive states refer to tropomyosin for which the accompanying actin binding sites are available for cross bridges to bind and generate force. Transitions between the states are governed by rate functions that depend on [Ca^{2+}]_{i} and *l*_{s}.

In the model of Courtemanche et al. (8), Ca^{2+} buffering by troponin is modeled by (19) where [Ca^{2+}]_{Trpn} is Ca^{2+}-bound troponin concentration, [Trpn]_{max} is total troponin concentration (70 μM) (8), and *K*_{m,Trpn} is half-saturation constant for troponin (0.5 μM) (8). In the model of Rice et al. (43, 44), the concentration of Ca^{2+} bound to high-affinity troponin sites is [HTRPNCa] and the dynamics are governed by (20) where [HTRPN]_{tot} represents the total troponin high-affinity site concentration and *k*_{htrpn}^{+} and *k*_{htrpn}^{−} are the Ca^{2+} on- and off-rates for troponin high-affinity sites (Table 1). The concentration of Ca^{2+} bound to low-affinity troponin sites is [LTRPNCa], and the dynamics are governed by (21) where [LTRPN]_{tot} represents the total troponin low-affinity site concentration and *k*_{ltrpn}^{+} and *k*_{ltrpn}^{−} are the Ca^{2+} on- and off-rates for troponin low-affinity sites (Table 1).

In our model, the Ca^{2+} transient is computed by the model of Courtemanche et al. (8) using an immediate formulation of Ca^{2+} binding by troponin (*Eq. 19*). The resulting Ca^{2+} transient is used to compute Ca^{2+} binding to troponin by *Eqs. 20* and *21*, and [LTRPNCa] is used to compute the tropomyosin rate from nonpermissive to permissive, as in the model of Rice et al. (43, 44). In the present study, we do not consider a feedback mechanism that influences the Ca^{2+} transient through a change in the affinity of troponin for Ca^{2+} binding as in *model 5* (43, 44). The choice between *model 4* and *model 5* is motivated in the discussion.

In *model 4*, the force generated by the sarcomeres depends on the fraction of tropomyosin in the force-generating states *N1*, *P1*, *P2*, and *P3*. We use the normalized force (F_{norm}), which is defined by (22) where *P1*_{max}, *P2*_{max}, and *P3*_{max} are defined as in Ref. 44 and φ(*l*_{s}) describes the physical overlap structure of thick and thin filaments within a sarcomere (44). When φ(*l*_{s}) = 1, all myosin heads are able to interact with actin in the single overlap zone; when φ(*l*_{s}) < 1, some of the filaments are in the double or nonoverlap zones. φ(*l*_{s}) is defined by (23)

In Fig. 2, the steady-state Ca^{2+}-force relation is presented for *model 4* (44). F_{norm} increases with increasing [Ca^{2+}]_{i} and with increasing *l*_{s}, with a maximum at *l*_{s} = 2.3 μm. To emphasize the dependence on [Ca^{2+}]_{i} and *l*_{s}, we will denote F_{norm} as a function: F_{norm}([Ca^{2+}]_{i},*l*_{s}).

### Mechanical Behavior of a Single Segment

The mechanical behavior of a single segment in our model is modeled as described by Solovyova et al. (49) by the classical three-element rheological scheme introduced by Hill in 1938 (14). Active force is generated by the contractile element (CE), and passive forces are generated in a series elastic element (SE) and a parallel elastic element (PE; Fig. 3). PE describes the force-length relation when the segment is not stimulated. CE and SE together describe the additional force generated on stimulation of the segment. The element lengths are *l*_{CE}, *l*_{SE}, and *l*_{PE}. The reference lengths, i.e., the lengths at which the segment is at rest and no force is applied, are *l*_{CE 0}, *l*_{SE 0}, and *l*_{PE 0}.

The force generated by the contractile element (F_{CE}) is defined by (24) where *f*_{CE} is a scaling factor, ν = −(d*l*_{s}/d*t*) represents the sarcomere shortening velocity, and F_{norm}([Ca^{2+}]_{i},*l*_{s}) is F_{norm} generated by the sarcomeres. The relation between the generated force and *v* is Hill's force-velocity relation (14, 17) and appears to be hyperbolic for skeletal and cardiac muscle (3, 9). We model the Hill relation by a function *f*_{ν}(ν) as proposed by Hunter et al. (17) (25) where ν_{max} is the maximum sarcomere shortening velocity and *c*_{ν} is a constant describing the shape of the hyperbole.

The forces generated in the SE and PE are nonlinearly dependent on their respective lengths *l*_{SE} and *l*_{PE} (49) and are defined by (26) and (27) where *l*_{SE 0} and *l*_{PE 0} denote the reference element lengths and *f*_{SE}, *k*_{SE}, *f*_{PE}, and *k*_{PE} are material constants describing the elasticity of the elements.

From mechanical equilibrium, it follows that F_{CE} must be equal to the force generated in the SE (F_{SE}). The total force generated by the segment (F_{segment}) is defined as F_{SE} + F_{PE}. Furthermore, *l*_{PE} must be equal to the *l*_{CE} + *l*_{SE} (Fig. 3). Therefore, during mechanical equilibrium (28) (29) (30) *l*_{CE}, *l*_{SE}, and *l*_{PE} are related to physiological sarcomere length (*l*_{s}) and reference sarcomere length (*l*_{s 0}) by *l*_{CE} = *l*_{s} and *l*_{CE 0} = *l*_{s 0} (49). The reference length of a segment is 0.01 cm and is related to *l*_{PE,0} by a scaling factor ξ. For segment *n*, we define the reference length *l*_{n 0} by (31) and the actual length *l*_{n} by (32) where *l*_{PE 0}^{n} and *l*_{PE}^{n} represent the reference length and the actual length of the PE of segment *n*. The λ for segment *n* (λ_{n}) is then defined by (33)

The parameters for the three-element mechanical model are obtained from Solovyova et al. (49) (Table 1). In Fig. 4, active force (F_{SE}), passive force (F_{PE}), and total force (F_{segment}) are presented for *l*_{s} = 1.7–2.5 μm and [Ca^{2+}]_{i} = 1.2, 0.9, 0.6, and 0.3 μM. When the sarcomeres generate force, i.e., *l*_{SE} > 0, *l*_{PE} = *l*_{CE} + *l*_{SE} is larger than *l*_{s} = *l*_{CE}. This results in a steeper increase of F_{PE} for increasing *l*_{s} and is in agreement with the passive force-length relation for intact cardiac muscle measured by Kentish et al. (23).

### Mechanical Behavior of a Cardiac Fiber

A cardiac fiber is modeled as a string of segments that are coupled in series. From mechanical equilibrium, it follows that the force F_{segment}^{n} generated by a single segment *n*, *n* ∈ *N*, is equal to the force generated by the fiber (F_{fiber}) (34) If we take into account that *l*_{n 0} may be different for each segment *n*, *n* ∈ *N*, the stretch ratio of the fiber (λ_{fiber}) is defined by (35) where *L* denotes the actual fiber length and *L*_{0} is the reference length.

In the present study, inhomogeneous cardiac tissue is represented by a 5-cm-long fiber with varying thickness and stiffness. The fiber is composed of 0.01-cm-long segments with 0.01- to 0.1-mm^{2} cross-sectional area. Tissue conductivity is related to stretch and may vary during the simulation. To enforce nonuniform stretch during the simulations, the diameter and stiffness of the left half of the fiber are varied, while the diameter (0.01 mm^{2}) and stiffness of the right half are normal. Linear interpolation is applied in a 0.5-cm transitional zone in the center of the fiber. Thick tissue is modeled by increasing the diameter of the segments, which affects the electrophysiological and mechanical properties of the tissue. Conductances, membrane surface, and the mechanical parameters *f*_{CE}, *f*_{SE}, and *f*_{PE} are scaled with the increase of the cross-sectional area. To simulate stiff tissue, the mechanical parameter *f*_{PE} is scaled. Scaling factors for maximum and minimum thickness are denoted by *t*_{max} and *t*_{min}, respectively, and scaling factors for maximum and minimum stiffness by *s*_{max} and *s*_{min}, respectively.

### Numerical Integration Scheme

To obtain criteria for the size of individual segments, we apply cable theory and consider subthreshold behavior along a fiber as previously described (30). For the bidomain parameters in Table 1, we obtain a length constant between 0.12 and 0.16 cm for λ = 1.0. When λ is increased to 1.4, the length constant decreases ∼15% for *R*_{junc 0}/*R*_{int 0} = 0.6. To obtain accurate simulation results, the fiber is modeled with segments that are 0.01 cm long, which is less than one-tenth of the length constant for λ ≤ 1.4. To solve the equations of the cellular bidomain model, we use a forward Euler scheme with a 0.01-ms time step to compute *V*_{mem} and an iterative method to solve the system of linear equations as described in Kuijpers et al. (30). Our method does not require matrix inversions and, therefore, is well suited to solve the system of equations when the conductivities change during the simulations as a result of stretch or contraction.

The ionic membrane currents are computed using a modified Euler method as described by Courtemanche et al. (8). To reduce computation time, the time step changes during the simulation as follows: a 0.01-ms time step is used shortly before and during the upstroke of the AP, and a 0.1-ms time step is used during repolarization and rest. The Ca^{2+}-force relation is computed using a forward Euler method with a fixed 0.1-ms time step, which is also the time step used to compute the cardiac mechanics (see appendix). Local conductivities are adjusted to the local λ whenever the mechanical state is updated.

To compare 0.1-ms (see above) with 0.01-ms time steps, we performed two simulations with the same parameter settings, but with different time steps. The differences in conduction velocity (θ), membrane currents, ionic concentrations, and mechanical forces were negligible, but computation time was reduced by 75%.

### Simulation Protocol

To illustrate the excitation-contraction coupling in our model, we performed single-cell simulations with constant *l*_{s} (isosarcometric contraction) and single-cell simulations with constant applied force (isotonic contraction). The influence of *I*_{sac} on the AP was investigated by application of a constant stretch to a single cell (isometric simulation). The cell was electrically stimulated with a frequency of 1 Hz. For investigation of spontaneous activity under stretch, simulations were performed with increasing stretch, but without electrical stimulation.

The influence of stretch on θ was investigated by stimulating the first segment of a 1-cm fiber. The fiber was short, such that contraction of early-activated segments did not affect impulse propagation in later-activated areas. The λ_{fiber} was kept constant during the simulation (isometric fiber contraction). We used longer (5 cm) fibers to investigate the influence of contraction on impulse propagation. Simulations were performed with contraction enabled and with contraction disabled. Disabled contraction was implemented by assuming that [Ca^{2+}]_{i} was equal to its resting value of 0.102 μM (8) when F_{norm}([Ca^{2+}]_{i},*l*_{s}) was computed. Thickness and stiffness were varied to simulate inhomogeneous cardiac tissue.

All simulations were performed over a 12-s period. Electrical stimulation was performed each 1 s (1 Hz) or each 0.5 s (2 Hz) by application of a stimulus current. In the case of single-cell simulations, a stimulus current of 20 pA/pF was applied for 2 ms as described by Courtemanche et al. (8). In the case of fiber simulations, the leftmost or the rightmost segment was stimulated by application of a stimulus current of 100 pA/pF until the membrane was depolarized. For the 1-cm fiber, the overall θ was measured by determination of the moment of excitation of two segments located 1 mm from each of the fiber ends. For the 5-cm (inhomogeneous) fiber, local θ was computed for each segment using the excitation time between two segments located 0.5 mm to the left and to the right in the nonstretched fiber.

## RESULTS

### Isosarcometric Contraction

Figure 5 illustrates the relation between the electrophysiology described by the model of Courtemanche et al. (8) and the force-producing states described by the model of Rice et al. (42–44). For the *V*_{mem} trace in Fig. 5, AP duration (APD) at 50% repolarization (APD_{50}) and APD at 90% repolarization (APD_{90}) are 184 and 304 ms, respectively. The AP amplitude and AP overshoot are 107 and 28 mV, respectively, and the maximum upstroke velocity [(d*V*_{mem}/d*t*)_{max}] is 187 V/s. For the Ca^{2+} transient, resting [Ca^{2+}]_{i} is 0.11 μM, peak [Ca^{2+}]_{i} is 0.87 μM, and time required to return [Ca^{2+}]_{i} to one-half of maximum [Ca^{2+}]_{i} is 178 ms. Since the dynamics of the concentration of Ca^{2+} bound to low-affinity troponin sites ([LTRPNCa]) are governed by a differential equation (21), the trace of [LTRPNCa] is smooth compared with that of [Ca^{2+}]_{i}.

From the traces of F_{norm} and individually normalized traces of F_{norm} for *l*_{s} = 1.7–2.3 μm in Fig. 6, it can be observed that peak force, time to peak force, and relaxation time increase with increasing *l*_{s}, which is consistent with the experimental data measured by Janssen and Hunter (18).

### Isotonic Contraction

In Fig. 7, traces of F_{norm}, F_{CE}, *l*_{s}, and *l*_{PE} are presented for simulations of isotonic contraction with applied force (F_{segment}) of 5–250 mN/mm^{2}. The AP and Ca^{2+} transient are the same as in Fig. 5. Less time is required to return F_{norm} to its resting value than in the case of isosarcometric contraction (Fig. 6, *top*). This is explained by shortening of the sarcomeres during contraction: a shorter sarcomere yields a lower contractile force (Fig. 2, *bottom*). The F_{CE} traces exhibit a plateau phase for F_{segment} ≤ 25 mN/mm^{2}. For F_{segment} = 250 mN/mm^{2}, *l*_{PE} remains constant, indicating no shortening.

### Effect of I_{sac} on AP

Figure 8 illustrates the effect of *I*_{sac} on the AP. *V*_{mem}, *I*_{sac}, *I*_{Ca,L}, and [Ca^{2+}]_{i} are presented for *I*_{sac} permeable to Ca^{2+} and *I*_{sac} not permeable to Ca^{2+} for λ = 1.00, 1.10, and 1.20. The cell was stimulated with a frequency of 1 Hz. With increasing λ, repolarization is prolonged and the resting *V*_{mem} is depolarized. *I*_{sac} is small during the plateau phase and larger during repolarization and rest, which is consistent with a reversal potential between 0 and −1 mV. *I*_{Ca,L} is somewhat lowered under stretch, and the Ca^{2+} transient is increased. The lowered *I*_{Ca,L} is explained by the Ca^{2+}-dependent inactivation of *I*_{Ca,L} (8). Interestingly, whether *I*_{sac} is permeable or not permeable to Ca^{2+}, the Ca^{2+} transient increases with increasing stretch. The characteristics for the AP and the Ca^{2+} transient are presented in Table 2 for *I*_{sac} permeable to Ca^{2+} and in Table 3 for *I*_{sac} not permeable to Ca^{2+}. For *I*_{sac} permeable to Ca^{2+}, peak [Ca^{2+}]_{i} and the time required for return of [Ca^{2+}]_{i} to one-half of maximum [Ca^{2+}]_{i} is increased ∼3%.

### Stretch-Induced APs

Figure 9 illustrates the effect of increasing λ in the presence of *I*_{sac}. *V*_{mem}, *I*_{sac}, *I*_{Ca,L}, and [Ca^{2+}]_{i} are presented for λ linearly increasing from 1.00 at 0-ms simulation to 1.25, 1.35, and 1.45 at 200-ms simulation; λ is constant after 200 ms. In both cases, stretch-induced APs are elicited for λ = 1.35 and 1.45. The APs for λ = 1.35 have a low upstroke steepness and are mainly driven by *I*_{Ca,L}. *I*_{Ca,L} increases faster for λ = 1.35 when *I*_{sac} is permeable to Ca^{2+}, which explains why *V*_{mem} reaches its maximum 50 ms earlier than when *I*_{sac} is not permeable to Ca^{2+}. For *I*_{sac} permeable to Ca^{2+} and for *I*_{sac} not permeable to Ca^{2+}, the sarcoplasmic Ca^{2+} flux signal for the Ca^{2+} release current (*I*_{rel}) is too small to trigger Ca^{2+} release from the SR. This explains why no Ca^{2+} transients are observed for λ = 1.35.

### Effect of R_{junc 0}/R_{int 0} on θ

To investigate the influence of *R*_{junc 0}/*R*_{int 0} on θ, we simulated impulse propagation along a 1-cm fiber for various *R*_{junc 0}/*R*_{int 0} and λ_{fiber}. *I*_{sac} was disabled in these simulations (*G*_{sac} = 0.0 μm/s). In Fig. 10, the overall θ is presented for *R*_{junc 0}/*R*_{int 0} = 0.0–1.0, and λ_{fiber} = 1.0–1.4. For *R*_{junc 0}/*R*_{int 0} = 1.0, θ is little affected by increasing λ_{fiber}; for lower values of *R*_{junc 0}/*R*_{int 0}, θ decreases with increasing λ_{fiber}. The decrease in θ is almost linear for *R*_{junc 0}/*R*_{int 0} = 0.4–0.8.

### Isometric Fiber Contraction

To investigate the influence of *I*_{sac} on impulse propagation, we simulated a series of isometric contractions in a 1-cm-long fiber. *G*_{sac} and λ_{fiber} were varied (*G*_{sac} = 0.0–0.020 μm/s, λ_{fiber} = 1.0–1.4). The leftmost segment was electrically stimulated with a frequency of 1 Hz. In Fig. 11, θ and (d*V*_{mem}/d*t*)_{max} are presented for various *G*_{sac} and λ_{fiber}. The decrease of θ with increasing *G*_{sac} is accompanied by a decrease of (d*V*_{mem}/d*t*)_{max}. When *I*_{sac} was permeable to Ca^{2+}, block of impulse propagation occurred for λ_{fiber} ≥ 1.35 for *G*_{sac} = 0.010 μm/s, λ_{fiber} ≥ 1.25 for *G*_{sac} = 0.015 μm/s, and λ_{fiber} ≥ 1.15 for *G*_{sac} = 0.020 μm/s. When *I*_{sac} was not permeable to Ca^{2+}, block of impulse propagation occurred for λ_{fiber} ≥ 1.25 for *G*_{sac} = 0.010 μm/s, λ_{fiber} ≥ 1.10 for *G*_{sac} = 0.015 μm/s, and λ_{fiber} ≥ 1.05 for *G*_{sac} = 0.020 μm/s.

In Fig. 12, the APs and traces of *I*_{sac}, *I*_{Ca,L}, and [Ca^{2+}]_{i} are presented for the center segment (λ_{fiber} = 1.00, 1.10, and 1.20). Traces of *I*_{sac,Na}, *I*_{sac,K}, and *I*_{sac,Ca} are shown in Fig. 13. As expected, *I*_{sac} increases with increasing λ during repolarization and rest. As in the single-cell simulations (Fig. 8), *I*_{Ca,L} decreases and [Ca^{2+}]_{i} increases with increasing λ. Since depolarization is mainly driven by *I*_{Na}, we further examine *I*_{Na} to investigate the cause of the decrease in θ and (d*V*_{mem}/d*t*)_{max}. The ionic current *I*_{Na} is defined by (36) where *G*_{Na} is the maximum *I*_{Na} conductance, *E*_{Na} is the equilibrium potential for Na^{+}, *m* is the fast activation variable, and *h* and *j* are the fast and slow inactivation variables (8). In Fig. 14, traces of *m*, *h*, *j*, and *I*_{Na} are presented for the center segment (λ_{fiber} = 1.00, 1.10, and 1.20). As λ_{fiber} increases, *h* and *j* are lower during rest and explain the lower *I*_{Na} during depolarization; i.e., the membrane is less excitable. Except for *I*_{to} and *I*_{Kur}, all ionic currents were similar during the upstroke and shortly after the upstroke (not shown). *I*_{to} and *I*_{Kur} were smaller and caused the less prominent notch of the AP (Fig. 12).

### Impulse Propagation Along a Homogeneous Fiber

To investigate the effect of contraction of early-activated areas on θ in later-activated areas, we simulated impulse propagation along a 5-cm-long fiber with λ_{fiber} = 1.00, 1.05, 1.10, 1.15, and 1.20 (*G*_{sac} = 0.015 μm/s, *P*_{Na}:*P*_{K}:*P*_{Ca} = 1:1:1). All simulations were performed with contraction enabled as well as with contraction disabled. Impulse propagation was initiated by application of a stimulus current to the leftmost segment.

In Fig. 15, traces of *V*_{mem}, *I*_{sac}, [Ca^{2+}]_{i}, and λ are presented for segments located 1.0, 2.5, and 4.0 cm from the stimulation site (λ_{fiber} = 1.15). With contraction enabled, the early-activated segments start contracting, so that λ increases for the later-activated segments, which results in an increased *I*_{sac} and a depolarized resting *V*_{mem}. In Fig. 16, θ and (d*V*_{mem}/d*t*)_{max} are presented with contraction disabled (*t*_{max} = 1, *s*_{max} = 1) and contraction enabled (*t*_{max} = 1, *s*_{max} = 1).

### Impulse Propagation Along an Inhomogeneous Fiber

To investigate the effect of inhomogeneity in tissue properties on θ, we simulated impulse propagation along an inhomogeneous 5-cm-long fiber with λ_{fiber} = 1.00, 1.05, 1.10, 1.15, and 1.20 (*G*_{sac} = 0.015 μm/s, *P*_{Na}:*P*_{K}:*P*_{Ca} = 1:1:1). For the left half of the fiber, *t*_{max} = 1–10 and *s*_{max} = 1–10. For the right half of the fiber, *t*_{min} = 1 and *s*_{min} = 1. Linear interpolation was applied in the central 0.5 cm of the fiber. As described above, all simulations were performed with contraction enabled as well as with contraction disabled. Impulse propagation was initiated by application of a stimulus current to the leftmost or the rightmost segment.

In Fig. 16, θ and (d*V*_{mem}/d*t*)_{max} are presented for various simulations after stimulation of the leftmost segment (λ_{fiber} = 1.15) with contraction disabled and with contraction enabled. In thick and/or stiff tissue (left half of the fiber), θ was larger; in the remaining, more stretched, parts, θ was smaller. In areas where the depolarization wave travels from thick tissue to thin tissue (*t*_{max} ≥ 5), was increased, which is explained by the smaller amount of charge required by the downstream segments to reach the excitation threshold. Decrease of (d*V*_{mem}/d*t*)_{max} and block of impulse propagation occurred in the inhomogeneous fibers when contraction was enabled.

In Fig. 17, θ and (d*V*_{mem}/d*t*)_{max} are presented for various simulations after stimulation of the rightmost segment (λ_{fiber} = 1.15) with contraction disabled and with contraction enabled. In thick and/or stiff tissue (left half of the fiber), θ was larger; in the remaining, more stretched, parts, θ was smaller. In areas where the depolarization wave travels from thin tissue to thick tissue (*t*_{max} ≥ 5), θ was decreased, which is explained by the larger amount of charge required by the downstream segments to reach the excitation threshold. Conduction block was not observed after stimulation from the right. Thus, for λ_{fiber} = 1.15, conduction block was unidirectional when contraction was enabled.

In Fig. 18, traces of *V*_{mem}, *I*_{sac}, [Ca^{2+}]_{i}, and λ are presented for three segments of an inhomogeneous fiber (λ_{fiber} = 1.15, *t*_{max} = 10, *s*_{max} = 1) with contraction enabled and with contraction disabled. As shown in Fig. 15, the early-activated segments start contracting, causing λ to increase for the later-activated segments, which leads to an increased resting *V*_{mem}. The AP of the segment at 4.0 cm had a low upstroke steepness and, similar to the AP for λ = 1.35 in Fig. 9, the Ca^{2+} transient was absent, such that no contraction occurred. From the rapid decrease in (d*V*_{mem}/d*t*)_{max} at ∼4.0 cm (Fig. 16, *bottom*), it can be concluded that this type of AP cannot generate enough current to propagate.

### Short Stimulation Intervals and Unidirectional Block

To investigate the effect of a shorter stimulation interval on impulse propagation, we stimulated the 5-cm fibers with an interval of 500 ms (2 Hz). In Table 4, θ for left stimulation at 1 Hz and at 2 Hz are presented for the left half (1.0 → 2.5 cm) and for the right half (2.5 → 4.0 cm) of the fiber. Since stimulation at 2 Hz can lead to alternating impulse propagation and conduction block, we distinguish between even (each 1 s) and odd (each 0.5 s) stimulation. The same data are presented in Table 5 for right stimulation.

From Tables 4 and 5, it can be concluded that stimulation at 2 Hz in general leads to slower conduction and conduction block at lower λ. This is explained by a longer ERP under stretch (Fig. 14). Figure 19 illustrates the subtle transition from conduction block in the leftmost 0.5 cm every other stimulation to normal impulse propagation every stimulation (λ_{fiber} = 1.05, *t*_{max} = 1, *s*_{max} = 1). In this case, the ERP decreased after each stimulation, such that after stimulation at 2,100 ms, the AP could propagate. The decrease in ERP is visible in Fig. 19 as the increasing *I*_{Na} inactivation gating variables *h* and *j* at the moment of stimulation (segment at 0.1 cm). After 2,100 ms, the cells in the fiber are stimulated at a higher frequency, which leads to a shorter APD and a more decreased ERP. Thus, impulse propagation at a higher frequency becomes a stable situation.

When an inhomogeneous fiber is stimulated from the right, conduction block may occur at lower λ. This can be explained by prolongation of the repolarization phase of the AP. Figure 20 illustrates this situation for λ_{fiber} = 1.10, *t*_{max} = 10, and *s*_{max} = 1. The extended repolarization phase of the segment at 4.0 cm (which is close to the stimulation site) is caused by contraction of the later-activated segments in the left half of the fiber.

## DISCUSSION

In our model, contraction of the cardiac fiber is triggered by the Ca^{2+} transient, which occurs after depolarization of the membrane. By modeling an *I*_{sac}, contraction of early-activated parts of the fiber leads to stretch in the later-activated parts and influences impulse propagation, APD, and ERP. For increasing levels of applied stretch, we observed conduction block, which can be unidirectional in an inhomogeneous fiber.

### Conduction Slowing and ERP

Our model provides two mechanisms to explain conduction slowing as observed in acutely dilated atria (6, 10, 11, 16, 40): *1*) the decrease in tissue conductivity due to stretch and *2*) a decreased membrane excitability caused by the *I*_{sac} (Fig. 11). In an experimental study, Eijsbouts et al. (11) reported a decreased θ and local conduction block when the right atrium of a rabbit was acutely dilated. They increased atrial pressure from 2 to 9 and 14 cmH_{2}O and measured λ as well as θ. With increasing pressure, θ first increases and then decreases for normal stimulation (240-ms interval). For fast stimulation (125-ms interval), θ decreases nonlinearly with increasing pressure (11). In our model, θ decreases linearly with increasing λ when no *I*_{sac} is present (Fig. 10). It is therefore likely that the nonlinear decrease in θ observed by Eijsbouts et al. is explained by a reduced excitability of the membrane, rather than a reduced tissue conductivity. Eijsbouts et al. (10, 11) also observed an increase in conduction block when the atrium was stimulated at a higher frequency, which is consistent with our findings (Tables 4 and 5).

Similar to our results, Shaw and Rudy (47) observed slowing of impulse propagation related to a reduced membrane excitability in a simulation study of impulse propagation in ischemic cardiac tissue. The extracellular K^{+} concentration ([K^{+}]_{o}) was increased, which leads to a depolarized *V*_{mem} and a reduced (d*V*_{mem}/d*t*)_{max}. Their simulation results (47) correspond to experimental results (19). Kléber and Rudy (27) explain the decreased θ in these experiments by a significant Na^{+} channel inactivation. The result is a depressed membrane excitability [reduced (d*V*_{mem}/d*t*)_{max}], reduced θ, and, eventually, conduction block (27). Thus their explanation for conduction slowing and block in ischemic tissue is similar to our explanation for conduction slowing and block under stretch.

In an experimental study, Sung et al. (51) observed a decrease in θ and an increase in APD when end-diastolic pressure was increased in the left ventricle of isolated rabbit hearts. Interestingly, the SAC blocker streptomycin had little effect on θ and APD (51). Satoh and Zipes (46) measured differences in ERP in the thin atrial free wall and the crista terminalis. Under stretch, the ERP of the thin atrial free wall was increased more than that in the thicker crista terminalis. Satoh and Zipes explain this difference by assuming that the thin free wall is more stretched than the thicker parts. Huang et al. (16) observed slow conduction related to a shorter stimulation interval in dilated atria, but they did not measure a significant change in atrial ERP after dilatation. Conduction slowing in our model at a 500-ms stimulation interval is attributed to a longer ERP under stretch. From these observations, we conclude that experimentally observed changes in conduction and atrial ERP can be explained by *I*_{sac}.

### Clinical Relevance

AF is associated with hemodynamic or cardiomechanical disorders such as hypertension, mitral valve disease, and cardiac failure (21). Ravelli et al. (41) found that atrial stretch caused by contraction of the ventricles influences atrial flutter cycle length in humans. Experimentally, it has been observed that acute atrial dilatation facilitates the induction and maintenance of AF in rabbit atria (2, 40) and in canine atria (16, 46). Bode et al. (2) report that the SAC blocker Gd^{3+} reduces the stretch-induced vulnerability to AF, confirming that *I*_{sac} plays a significant role in the vulnerability to AF in acutely dilated atria.

In the present study, we observed conduction slowing, an increased ERP, and (unidirectional) conduction block with increasing stretch. These phenomena are attributed to *I*_{sac} and can lead to alternating impulse propagation and contractions at a stimulation frequency of 2 Hz. Conduction slowing (16), unidirectional block (27), and dispersion in atrial ERP (46) are related to the inducibility of AF. In the present study, these effects begin at λ = 1.15 and 1.05 for stimulation at 1 and 2 Hz, respectively. Bode et al. (2) gradually increased intra-atrial pressure in rabbit hearts up to 30 cmH_{2}O. They could not induce AF in the undilated atrium at 0 cmH_{2}O, but they observed a 50% probability of AF induction at 8.8 cmH_{2}O (baseline) and at 19.0 cmH_{2}O (after Gd^{3+}), which increased to 100% (baseline) and 90% (after Gd^{3+}) when pressure was further increased (2). However, stretch was not measured in that study. Eijsbouts et al. (11) measured λ = 1.16 ± 0.14 at 9 cmH_{2}O, confirming our model predictions that the vulnerability to AF is substantially increased when λ is ∼1.15.

### Model Validity and Limitations

To our best knowledge, our model is the first to integrate cardiac electrophysiology and cardiomechanics with physiological details such as ionic membrane currents, intracellular Ca^{2+} handling, and cross-bridge formation. In our model, changes in impulse propagation under stretch are related to *I*_{sac} and to a reduced conductivity. We do not consider other mechanisms that could influence impulse propagation, such as stretch-related function of other membrane channels, autonomic reflexes, and metabolic changes.

The validity of our model largely depends on the validity of the underlying models and parameters. Validity and limitations of the models for the ionic membrane currents, cross-bridge formation, and cardiomechanics are extensively discussed elsewhere (8, 44, 49, 58). Here, we discuss the validity and limitations of the integrated model with respect to the Ca^{2+}-force relation, excitation-contraction coupling, stretch and fiber conductivity, *I*_{sac}, intracellular ion concentrations, Ca^{2+}-troponin binding, and the Ca^{2+} transient.

#### Ca^{2+}-force relation.

Rice et al. (44) proposed five models of isometric force generation in cardiac myofilaments. These are constructed assuming different subsets of three putative cooperative mechanisms (44). *Model 4* assumes that the binding of a cross bridge increases the rate of formation of neighboring cross bridges and that multiple cross bridges can maintain activation of the thin filament in the absence of Ca^{2+}. The model also simulates end-to-end interactions between adjacent troponin and tropomyosin. The hypothesis that cross-bridge binding increases the affinity of troponin for Ca^{2+} is assumed by *model 5*, but not by *model 4*.

To choose between *model 4* and *model 5* for the present study, we studied the Ca^{2+}-force relation (see Fig. 2 for *model 4*) and the isosarcometric twitches (see Fig. 6 for *model 4*). We found better agreement between the Ca^{2+}-force relation obtained by *model 4* and the experimental results of Kentish et al. (23), in particular for >1.9-μm-long sarcomeres. When we compared the isosarcometric twitches, we found that, for *model 5*, the peak force was lower and the latency to peak force was increased for longer sarcomeres. Compared with the experimental data measured by Janssen and Hunter (18), the latency to peak force increased too greatly with sarcomere length. Our findings confirm the finding by Rice et al. (44) that the hypothesis that cross-bridge binding increases the affinity of troponin for Ca^{2+} is not crucial to reproduce the experimental results. Since the twitches obtained by *model 4* better resemble the experimental results from Janssen and Hunter, we have chosen *model 4* to describe the Ca^{2+}-force relation.

#### Excitation-contraction coupling.

To compute F_{norm} during contraction, Rice et al. (44) used a Ca^{2+} transient with a peak [Ca^{2+}]_{i} of 0.97 μM and 130 ms to return [Ca^{2+}]_{i} to half of maximum [Ca^{2+}]_{i}. Our peak forces are smaller and the traces are less prolonged than those of Rice et al. (44). This is explained by the lower peak value of [Ca^{2+}]_{i} obtained from the model of Courtemanche et al. (8). However, the main characteristics, i.e., increasing peak force, increasing time to peak force, and increasing relaxation time with increasing *l*_{s}, are observed in our model and are in agreement with experimental measurements (18). These characteristics are important with respect to the Frank-Starling mechanism, which states that when the amount of blood flowing into the heart increases, the wall becomes more stretched and the cardiac muscle contracts with increased force. Furthermore, isosarcometric twitch duration (Fig. 6, *top*) is longer than isotonic twitch duration (Fig. 7, *top*). This is explained by shortening of the sarcomeres during isotonic contraction and is in agreement with experimental results (17).

#### Stretch and fiber conductivity.

In our model, *g*_{int} and *g*_{ext} are determined by λ on the basis of the assumption that 60% of *R*_{int} is attributed to *R*_{junc}. The assumption that *R*_{int} is determined by *R*_{myo} only is equivalent to *R*_{junc 0}/*R*_{int 0} = 0.0 and would lead to a faster decrease in θ for increasing λ (Fig. 10). Since intercellular coupling in cardiac tissue is through the gap junctions and since the number of gap junctions does not change when the tissue is acutely stretched, we believe that neglecting the distinction between *R*_{myo} and *R*_{junc} leads to overestimation of the effect of λ on conductivity.

#### I_{sac}.

We model *I*_{sac} as a nonselective cation current with *E*_{sac} = 0 to −1 mV. Our model of *I*_{sac} has a near-linear current-voltage relation, which can be approximated by (37) with *g*_{sac} = 0.0027 nS/pF for λ = 1.0, *g*_{sac} = 0.0049 nS/pF for λ = 1.2, and *g*_{sac} = 0.0088 nS/pF for λ = 1.4. Wagner et al. (57) used an externally applied current with *g*_{sac} = 0.0083 nS/pF and *E*_{sac} −10 mV to elicit an AP in atrial cells from the rat. In our model, stretch-induced APs are elicited for λ = 1.35 (Fig. 8), which corresponds to *g*_{sac} = 0.0076 nS/pF and is in a similar range.

#### Intracellular ion concentrations.

The ionic currents of the model of Courtemanche et al. (8) interact with [Na^{+}]_{i}, [K^{+}]_{i}, and [Ca^{2+}]_{i} (8). To model the influence of *I*_{sac} on the intracellular ion concentrations, we assume that SACs are permeable to Na^{+}, K^{+}, and Ca^{2+}. Whether SACs in human atrial cells are permeable to Ca^{2+} is still a matter of debate. To investigate the effect of permeability of SACs for Ca^{2+}, we performed simulations in which *I*_{sac} was permeable to Ca^{2+} and not permeable to Ca^{2+}. Our results show an increase in [Ca^{2+}]_{i} under stretch, whether *I*_{sac} is permeable to Ca^{2+} or not permeable to Ca^{2+} (Fig. 8). However, when *I*_{sac} is permeable to Ca^{2+}, [Ca^{2+}]_{i} is increased ∼3% compared with *I*_{sac} not permeable to Ca^{2+} (Tables 2 and 3). Our findings are in agreement with experimental observations that [Ca^{2+}]_{i} increases in response to stretch (4, 52). Kamkin et al. (20) observed that the behavior of SACs in isolated human atrial cells did not change when a Ca^{2+}-free external solution was used, suggesting that the current through the SACs was preferentially carried by Na^{+}, rather than by Ca^{2+} (57). In Fig. 13 (*left*), the current density of *I*_{sac,Ca} is only 5% of the current density of *I*_{sac,Na}, which explains the finding by Kamkin et al.

#### Ca^{2+}-troponin binding.

Binding of Ca^{2+} to troponin is modeled in two different ways: in the model of Courtemanche et al. (8), an immediate formulation (*Eq. 19*) to describe Ca^{2+} binding to troponin is used; in the model of Rice et al. (44), binding of Ca^{2+} to troponin is modeled by *Eqs. 20* and *21*. The resulting [LTRPNCa] is then used to compute the force generated by the sarcomeres. Although having two different ways to model binding of Ca^{2+} to troponin is not elegant, one cannot simply replace the immediate formulation proposed by Courtemanche et al. with the formulation proposed by Rice et al. This would be a major modification of the model of Courtemanche et al. requiring adaptation of various parameters related to the Ca^{2+} fluxes and ionic currents, which is beyond the scope of the present study.

#### Ca^{2+} transient.

As shown in Figs. 9 and 18 (*left*), although there is an AP, the Ca^{2+} transient is absent. The APs for which no Ca^{2+} transient occurred are characterized by a low steepness of the AP upstroke [low (d*V*_{mem}/d*t*)_{max}] and were mainly driven by *I*_{Ca,L}. The absence of the Ca^{2+} transient is explained by the absence of Ca^{2+} release from the junctional SR, which may be caused by the lower (d*V*_{mem}/d*t*)_{max}. In our model, the effect of the intercellular currents on the ion concentrations is not taken into account. In case the AP is mainly driven by *I*_{Ca,L}, [Ca^{2+}]_{i} in the downstream cell may slowly rise because of the Ca^{2+} in the current flow between the cells. This increase in [Ca^{2+}]_{i} may then trigger buffered Ca^{2+} release. One may expect that a similar effect can be obtained when the SACs are permeable to Ca^{2+}. Although AP rises faster in the presence of *I*_{sac,Ca} (λ = 1.35 in Fig. 9), a Ca^{2+} transient was not observed. In simulations where Ca^{2+} transients did occur, the effect of *I*_{sac,Ca} was limited to an ∼3% increase in [Ca^{2+}]_{i} (Tables 2 and 3). On the basis of these observations, we conclude that the permeability of *I*_{sac} for Ca^{2+} contributes little to the changes in electrophysiological behavior under stretch.

In conclusion, in our model, conduction slowing and block are related to the amount of stretch and are enhanced by contraction of early-activated segments. Conduction block can be unidirectional in an inhomogeneous fiber and is promoted by a shorter stimulation interval. Our observations are in agreement with experimental results and provide an explanation for the increased inducibility of AF observed in acutely dilated atria.

## APPENDIX

#### Solving the three-element model.

The numerical scheme to solve the equations for the three-element mechanical model is based on the scheme described by Solovyova et al. (49). The forces and lengths of CE, SE, and PE are computed by introducing *l*_{1} = *l*_{CE} − *l*_{CE 0} and *l*_{2} = *l*_{PE} − *l*_{PE 0}. Furthermore, it is assumed that *l*_{PE 0} = *l*_{CE 0}, from which it follows that *l*_{SE 0} = *l*_{PE 0} − *l*_{CE 0} = 0 μm. *Equations 26* and *27* can now be rewritten as (38) and (39) The mechanical state of each segment is now defined by *l*_{1}, *l*_{2}, d*l*_{1}/d*t*, d*l*_{2}/d*t*, and F_{norm}([Ca^{2+}]_{i},*l*_{s}). Each simulation time step F_{norm} is computed using [Ca^{2+}]_{i} obtained from the model of Courtemanche et al. (8) and *l*_{s} of the former time step. Next, *l*_{1} and *l*_{2} are updated using a forward Euler step and d*l*_{1}/d*t* and d*l*_{2}/d*t* of the former time step. New values of d*l*_{1}/d*t* and d*l*_{2}/d*t* are then computed using *Eqs. 24* and *25* (40) from which sarcomere shortening velocity (ν) can be obtained by (41) F_{CE} can be obtained from *Eqs. 38* and *28* by (42) Since *l*_{s} = *l*_{CE} and *l*_{1} = *l*_{CE} − *l*_{CE 0}, we obtain for ν (43) from which d*l*_{1}/d*t* follows immediately.

For isometric single-segment simulations, d*l*_{2}/d*t* = 0 and the generated force can be directly computed from *Eqs. 38* and *39*. For isotonic simulations, d*l*_{2}/d*t* can be obtained from dF_{segment}/d*t* using F_{segment} = F_{SE} + F_{PE} and *Eqs. 38* and *39* (44) and, by taking the derivative (45) from which d*l*_{2}/d*t* can be computed by (46) Note that during isotonic contraction dF_{segment}/d*t* = 0.

#### Computing cardiac fiber mechanics.

To obtain a solution for a multiple-segment simulation, we define α and β by (47) and (48) *Equation 46* can be formulated for each segment *n* by (49) where α_{n} and β_{n} denote α and β for segment *n*. Using *l*_{2}^{n} = *l*_{PE}^{n} − *l*_{PE 0}^{n} and *l*_{n} = ξ_{n}*l*_{PE}^{n} (*Eq. 32*), we obtain for fiber length *L* (50) Using F_{segment}^{n} = F_{fiber} for each segment *n* ∈ *N* (*Eq. 34*) and introducing *a* and *b*, we obtain (51) where *a* and *b* are defined by (52) and (53) Finally, using the definition of λ_{fiber} = *L*/*L*_{0} (*Eq. 35*) and *Eq. 51*, we obtain (54) from which follows (55) For isotonic simulations, F_{fiber} is constant and dF_{fiber}/d*t* = 0. Using F_{segment}^{n} = F_{fiber} (*Eq. 34*), d*l*_{2}^{n}/d*t* can be obtained from *Eq. 46* for each segment *n* ∈ *N*. For isometric simulations, λ_{fiber} is constant and dλ_{fiber}/d*t* = 0. In the isometric case, F_{fiber} and its derivative dF_{fiber}/d*t* can be obtained from *Eq. 55*.

In summary, the mechanical state of segment *n* is described by *l*_{1}^{n}, *l*_{2}^{n}, d*l*_{1}^{n}/d*t*, d*l*_{2}^{n}/d*t*, and F_{norm}^{n}. F_{norm}^{n} is dependent on [Ca^{2+}]_{i} obtained from the model of Courtemanche et al. (8) and the sarcomere length *l*_{s} = *l*_{s 0} + *l*_{1}^{n}. For each time step, *l*_{1}^{n} and *l*_{2}^{n} are computed using a forward Euler step, and d*l*_{1}^{n}/d*t* and d*l*_{2}^{n}/d*t*, respectively. However, initial values for *l*_{1}^{n} and *l*_{2}^{n} remain to be defined. Initially, it is assumed that the electrophysiological state of all segments is resting; i.e., *V*_{mem} is −81 mV and [Ca^{2+}]_{i} is 0.102 μM (8). For such low [Ca^{2+}]_{i}, F_{norm} is small and we assume F_{norm}^{n} = 0 for each segment *n*. Thus, F_{CE}^{n} = 0, and since F_{SE}^{n} = F_{CE}^{n}, the force generated by the segment must come from the PE, i.e., F_{PE}^{n} = F_{segment}^{n}. Since F_{SE}^{n} = F_{CE}^{n} = 0, *l*_{1}^{n} = *l*_{2}^{n} (*Eq. 38*). For homogeneous tissue, the material properties *k*_{PE} and *f*_{PE} are equal for all segments. For isotonic simulations, *l*_{2}^{n} can be obtained from *Eq. 39* by (56) For isometric simulations, *l*_{2}^{n} can be directly obtained from the initial stretch ratio λ_{fiber} by (57) For inhomogeneous tissue, it is assumed that, initially, λ_{fiber} = 1 and F_{fiber} = 0. In that case, *l*_{1}^{n} = *l*_{2}^{n} = 0 for all segments *n*. During the simulation, stretch (isometric simulation) or force (isotonic simulation) is slowly increased until the desired values are reached. This process typically requires 200 ms of simulation. Finally, it is assumed that, initially (58) for homogeneous and inhomogeneous tissue.

## Footnotes

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 © 2007 by the American Physiological Society