Arrhythmogenic Current Generation by Myofilament-Triggered Ca2+ Release and Sarcomere Heterogeneity.

Heterogeneous mechanical dyskinesis has been implicated in many arrhythmogenic phenotypes. Strain-dependent perturbations to cardiomyocyte electrophysiology may contribute to this arrhythmogenesis through processes referred to as mechanoelectric feedback. Although the role of stretch-activated ion currents has been investigated using computational models, experimental studies suggest that mechanical strain may also promote arrhythmia by facilitating calcium wave propagation. To investigate whether strain-dependent changes in calcium affinity to the myofilament may promote arrhythmogenic intracellular calcium waves, we modified a mathematical model of rabbit excitation-contraction coupling coupled to a model of myofilament activation and force development. In a one-dimensional compartmental analysis, we bidirectionally coupled 50 sarcomere models in series to model calcium diffusion and stress transfer between adjacent sarcomeres. These considerations enabled the model to capture 1) the effects of mechanical feedback on calcium homeostasis at the sarcomeric level and 2) the combined effects of mechanical and calcium heterogeneities at the cellular level. The results suggest that in conditions of calcium overload, the vulnerable window of stretch-release to trigger suprathreshold delayed afterdepolarizations can be affected by heterogeneity in sarcomere length. Furthermore, stretch and sarcomere heterogeneity may modulate the susceptibility threshold for delayed afterdepolarizations and the aftercontraction wave propagation velocity.

Significance

Cardiac arrhythmia is a major contributor to mortality and can result from spatiotemporal perturbations to the electrical, mechanical, or biochemical properties of the myocyte. Nonhomogeneous wall stresses give rise to nonuniform tissue strains and sarcomere lengths that are thought to contribute to these arrhythmogenic alterations and, in certain cases, to be directly involved in the initiation of arrhythmia via mechanoelectric feedback. The attention given to mechanoelectric feedback mechanisms has largely been concentrated on stretch-activated currents. However, mechanical perturbations to myofilament calcium buffering are also considered to be of importance. Using a strongly coupled electromechanical myocyte model, this study suggests that heterogeneous sarcomere lengths and sarcomere strains may contribute to the initiation of propagating calcium waves and modification of their propagation velocity.

Introduction

Cardiac arrhythmias are often initiated by spatiotemporal instabilities in electrical activity that can be triggered by various mechanisms including localized mechanical disturbances. The stability and maintenance of arrhythmias is affected by cardiac conduction properties that can also be modified by changes in myocyte mechanical loading and strain (1,2).

Elementary calcium release events from the sarcoplasmic reticulum (SR), called calcium sparks, can result in calcium transients (CaTs), but only when triggered spatially uniformly and simultaneously during normal action potentials (APs) and under calcium-overloaded conditions (3,4). Under normal conditions, calcium sparks do not propagate either spontaneously (5, 6, 7) or evoked (8, 9, 10, 11, 12), whereas in calcium-overloaded cells, the sites of the calcium sparks have been linked to the origin of calcium waves (3,5,13). Increased calcium spark frequency together with simultaneous recruitment of these spark-initiated calcium waves under unstimulated conditions is thought to be the origin of oscillatory membrane potentials that contribute to triggered arrhythmias (4,14, 15, 16, 17).

Although the electrical properties and origin of spark-induced calcium waves are well studied (3,5,13,18, 19, 20) and have been associated with asynchronous myofilament motion (21), the role of abnormal mechanical events and their contributions to modifying calcium waves are less well known. Experiments suggest that transient mechanical unloading can induce calcium release from myofilaments by dissociating calcium from troponin C (TnC) (22, 23, 24, 25, 26, 27). In disease, the inhomogeneous mechanical properties of the myocardium, together with abnormal SR calcium handling, are associated with increased incidence of calcium waves that may cause not only aftercontractions but also transient inward currents and arrhythmogenic delayed afterdepolarizations (DADs) (1,28). The experimental data reported by Wakayama et al. (22) of quick stretch and release experiments on cardiac trabeculae suggest that the magnitude of calcium waves increases at greater stretch-releases as a result of calcium release from TnC during rapid unloading. In turn, the locally induced increase in intracellular calcium triggers calcium release from the adjacent SR as the initiating trigger for calcium waves. During the calcium wave, the amount of calcium released from the SR is dictated by the magnitude of the initial trigger. Therefore, the local calcium release from the SR determines the calcium amount released in the adjacent sites, modifying the amplitude of the resulting DAD. Thus, mechanically induced arrhythmogenic calcium waves require not only a region of calcium-loaded SR for the propagation of the calcium wave but also a mechanical event that is needed as the trigger of the wave. ter Keurs et al. (23) found that local heterogeneity of sarcomere shortening can be the origin of calcium waves in experiments in which they induced a highly localized discontinuity in sarcomere shortening by subjecting isolated trabeculae to a fine laminar jet flow of butane-dione monoxime in an electromechanical decoupler. Most importantly, fast-propagating calcium waves have also been observed when stretch-activated currents (SACs) were blocked and extracellular calcium concentration was increased, as reported in several studies (22,23,26,29,30).

Previous computational models investigated the triggers of myofilament-induced calcium waves by simulating the effects of stretch-release as strain-rate-dependent changes in calcium sensitivity of the myofilaments, the role of SACs, and calcium release events as mediators of mechanically stimulated electrical instabilities. In contrast, this study investigates the propagation of mechanically triggered waves (23,31). Specifically, we assessed how specific calcium-myofilament interactions affect the induction of calcium waves and their propagation velocity.

Materials and Methods

Electromechanical model

We developed a one-dimensional (1D) strongly coupled computational electromechanics (EM) model of the rabbit ventricular cardiomyocyte to investigate possible mechanisms of triggered propagating contractions (TPCs) and the determinants of their wave speed.

For the electrophysiology (EP) model, we used an implementation of the established rabbit-specific AP ventricular myocyte model by Shannon et al. (32) (distributed by University of California Davis; https://somapp.ucdmc.ucdavis.edu/Pharmacology/bers/). The model represents calcium-induced calcium release (CICR) from the SR at the dyadic cleft to reflect the calcium gradients generated during normal excitation-contraction coupling and how those gradients interact with the localization of calcium-sensitive transporters (33) as shown in Fig. 1. In this model, although calcium can bind to both high- and low-affinity binding sites on TnC, magnesium competes with calcium only for the high-affinity site (32).

Figure 1

Schematic of cardiac excitation-contraction coupling as presented by Bers (70), showing the calcium transport in a ventricular myocyte with the sarcolemma, subsarcolemmal compartment (SS), the L-type calcium channel (L-type), the sarco/endoplasmic reticulum calcium-ATPase (SERCA), and the sarcoplasmic reticulum (SR). The red arrows indicate calcium influx resulting in an increase in cytosolic calcium concentration, leading to contraction of the myofilaments, which is the feedback mechanism of mechanoelectric feedback (MEF) from EP to mechanics. The black arrows show the pathway for calcium removal from the cell. In particular, the black arrow from the myofilaments shows the calcium release from TnC into the cytosol during relaxation, which is the feedforward mechanism of MEF from mechanics to EP. The time course of an AP (blue line), a CaT (red line), and contraction (force) (green line) are shown at the bottom. RyR, ryanodine receptor; TnC, troponin C.

For myofilament mechanics, the model by Rice et al. (34) for rabbit cardiac muscle contraction was used (downloaded from CellML; https://models.cellml.org/exposure/a31519a27f4c2fe6158e04fd40eeda98/rice_wang_bers_detombe_2008.cellml/cellml_codegen). It describes the activation of the thin filament by intracellular calcium binding to TnC and cross-bridge (XB) cycling. To simulate stretch-dependent calcium sensitivity of the myofilaments, the binding of calcium to TnC relies on mechanical perturbations imposed by differences in sarcomere resting length and changes in cell loading. These mechanical influences only impact the occupancy of the low-affinity TnC states in the EP model, with strong coupling of EP and sarcomere mechanics.

Considering that ryanodine receptor (RyR) clusters and junctions are spaced ∼2 μm apart longitudinally over an ∼100 μm long cardiomyocyte (35), to obtain a 1D representation of a myocyte, 50 EM models were coupled in series by calcium fluxes in the longitudinal direction, as well as the mechanical interaction of adjacent sarcomeres. The fractional volume of each subcellular compartment was chosen as reported in Shannon et al. (32), assuming that the myocyte is separated into 50 uniform, equidistant units (Nunit = 50). Therefore, the volume of each unit was assumed to be th of the volume of the myocyte, i.e., Vunit = 2 × π × 10.252 × 10 μm3 ×  = 4 × π × 10.252 μm3.

Incorporation of TPC triggers

To empirically approximate the triggers of TPCs as strain-dependent changes in calcium affinity to the myofilaments, we parameterized a zero-dimensional (0D) representation of the (mechanical) model to reproduce the CaT and twitch dynamics as reported by Wakayama et al. (22) (see Fig. 2 B, CaT; Fig. 2 D, force; and Fig. 2 F, stretch protocol). In particular, we used the parameter-fitting algorithm proposed by Bueno-Orovio et al. (36) to adjust the on and off rates for calcium binding to TnC (kon and koff, respectively; see Eq. 1) and to alter the strain-dependence of the XB cycling in the Rice model (34). During the same optimization, all four constants were reformulated to have an exponential dependence on strain-rate (). Similarly to our previous study (31), we fitted the 0D representation of the model to experimental data (22) by minimizing the mean-square error (MSE) in magnitudes for the twitch force and calcium wave in TPCs.

Figure 2

Simulation data and the corresponding experimental data of Wakayama et al. (22) (top) and Takamatsu et al. (19) (bottom), which were employed for the model parameterization using the parameter-fitting algorithm of Bueno-Orovio et al. (36). (Top) The effect of square waveform stretch experiments at 10% stretch (red line), 5% stretch (green line), control (dashed black line), and 10% compression (blue line) are shown on the (A) simulated CaT, (B) experimentally measured CaT, (C) simulated normalized force (here, normalized force refers to the maximal normalized force given full activation as described in (34)), and (D) experimentally measured force. (E) Simulated square waveform stretch experiments and (F) experimentally recorded square waveform stretch experiments are shown. The solid black line in (A) shows the CaT in the original, uncoupled Shannon et al. model (32) and in (C) the normalized force for the original, uncoupled Rice et al. model (34). (G) Similar to the experimental data reported in Fig. 2 C of (19) (dashed blue lines), the computational model (sarcomere resting length of 1.89 μm) reproduces a linearly propagating calcium wave through the myocyte with a velocity of ∼100 (solid red lines). (H) The calcium wave front at spatially different locations within the longitudinal diameter of a cardiac myocyte is shown for experiments reported in Fig. 2 D of (19) (dashed blue line) and simulations (solid red line). To see this figure in color, go online.

The algorithm-based changes resulted in a reduction of the scaling factors σn and σp, which influence the strain dependency of the transition rate from the post-force-generating state to the permissive state in the XB cycle. As a result of the strain dependency of the transition rate, not only is more calcium bound and then released by the myofilaments during rapid stretches, but the force development is also increased. Hence, the impact of stretch on the XB cycling rate is enhanced in the computational model, showing similar cytosolic calcium fluctuations induced by acute square waveform stretches as reported in experiments (22), shown in Fig. 2 (Fig. 2 A, simulated CaT; Fig. 2 C, simulated force; and Fig. 2 E, simulated stretch). Testing the main characteristics of the model as presented in Rice et al. (34) (steady-state force-sarcomere length relations, steady-state force-calcium relations (force-calcium relations) including sarcomere length effects, steady-state sarcomere length-calcium relations for unloaded cells, steady-state force-velocity relations, isometric twitches including calcium activation and sarcomere length effects, cell shortening twitches as a function of activator calcium, and effects of sarcomere length control on the intracellular CaTs; data not shown), the differences between the original and the updated mechanics model were <1% when comparing the MSE. The resulting updated formulation for the strain-dependent changes in calcium affinity to the myofilaments can be formulated as

These equations are used for the calcium buffering on TnC, in particular TnC (low-affinity sites) (TnCL) and TnC (high-affinity sites) (TnCH).with Cacyto the cytosolic calcium concentration, heaviside(v) the Heaviside step function of MATLAB (The MathWorks, Natick, MA), x0 = 0.007 μm strain induced by head rotation, xXBPostR the mean distortion of the post-force-generating state, and the transition rate from the post-force-generating state to the permissive state of the three-state Markov model of the XB cycle.

Square waveform stretch experiments shown in Fig. 2 E with 5 and 10% stretch (green and red lines, respectively) result in an increase of the peak CaT during stretch-release by 22 and 53% from control (black lines), respectively (Fig. 2 A), compared to 11 and 32% in the experimental data reported in Fig. 1 of Wakayama et al. (22) and illustrated in Fig. 2 B. Similarly, the peak force increased in the simulations by 24 and 47% from control, respectively (Fig. 2 C), compared to 8 and 27% as reported by Wakayama et al. (22) (Fig. 2 D). In contrast to stretch, compression square waveform experiments with 10% compression decreased the force peak by 46% from control and the myofilament-triggered CaT by 9% (Fig. 2, A and C, blue lines) compared to 7 and 20%, respectively, reported in the experimental data (Fig. 2, B and D).

Time-dependent calcium diffusion

For every unit i, calcium can diffuse from its left neighbor (i − 1) and right neighbor (i + 1) as proposed by Voigt et al. (37) (see Fig. 3, black arrows). Calcium diffusion was implemented for three compartments: SR, subsarcolemma (SS), and cytosol (albeit with different diffusion time constants). For intraunit diffusion in the longitudinal direction in the 1D model, we adapted the existing 0D model formulation for membrane ion kinetics from Shannon et al. (32) by adjusting the different diffusion rates by a factor of change of area between the compartments due to segmentation. The diffusion time constants governing calcium flux within and between the Nunit sarcomere units were adjusted to 1) retain the calcium transient magnitudes in each sarcomere compartment as published in the original ionic model (32) and 2) reproduce calcium wave velocities (100 ).

Figure 3

Schematic for two of the 50 units of the 1D model. The black arrows indicate the calcium diffusion between the SS compartments (top) described by the diffusion time constant τSS, the SR compartments (middle) described by the diffusion time constant τSR, and the cytosolic compartments (bottom) described by the diffusion time constant τSR. The red arrow shows how changes in length of the sarcomere in unit i (δSLi) affects the length of the adjacent sarcomere in unit (i + 1) (δSLi). RyR, ryanodine receptor.

Some experimental studies have reported a wave velocity of 100 at extracellular calcium concentrations (Cao) of 1.8–2 mM (19,38) (see Fig. 2, G and H, dashed blue lines), whereas others report velocities of 72 at Cao = 2 mM, 80 at Cao = 5 mM, and 90 at Cao = 15 mM. In our model, propagating calcium waves only occurred at Cao ≥ 4 mM, with the final parameterized model producing a calcium wave velocity of 100 at Cao = 4 mM (see Fig. 2 G). In particular, we used the parameter-fitting algorithm of Bueno-Orovio et al. (36) to adjust the parameterization for the diffusion time constants reported in (37) so that it replicates a reasonable gain of calcium-mediated DADs and appropriate intracellular propagation of the calcium waves (3,39) in ventricular cardiomyocytes. Locally in each unit, the CaTs have a peak magnitude of 1.6 μm, with a slight downward trend that results from the slowing of the calcium diffusion from the unit in which the wave was initiated to the boundaries. Similarly, the experimental data from Takamatsu et al. (19) suggest that the calcium waves rise steeply at their front (illustrated in Fig. 2 H, blue dashed line), similar to the CaT in our model as shown in Fig. 2 H (red line). To initiate a calcium wave, the open probability of the RyRs in the first unit of the myocyte was increased to 100% for 1 ms.

For all compartments except on the boundary, we have the net diffusion of calcium between a unit i and its adjacent units:

For the boundary units, we have the following boundary conditions:with τcyto = 1.2 ms, τSR = 150 ms, and τSS = 0.4533 ms.

Thus, the calcium concentrations for the different compartments within each unit can be updated as follows (see Shannon et al. (32) for detailed description):with , k ∈ [SS, cyto, SR], the calculation for the calcium concentrations as in Shannon et al. (32).

Owing to the almost instantaneous propagation of the electrical signal, the membrane potential was calculated using the average of the total membrane current in each unit,with Ii the total membrane current in unit i ∈ [1, …, 50] and Cm the membrane capacitance.

Incorporation of changes in interunit distances

For the sarcomeric model, we used the mean-field three-state Markov model approach of Rice et al. (34) (rabbit representation) to represent XB cycling, which represents biophysical mechanisms of the myofilaments with ordinary differential equations (ODEs). The primary cooperative mechanism in the Rice model is the XB-release unit cooperativity to generate the steep force-calcium relationship.

In our spatially explicit model, one sarcomere model based on (34) was allocated to each EP unit. To strongly couple neighboring myofilament contractile kinetics, adjacent sarcomeres affected each other by their fractional change in length (see Fig. 3, red arrow):with the internal, isotonic change (as calculated in (34)) of the ith SL for i ∈ {2, ..., 49} calculated independently from its neighbors by using the regional calcium concentration.

At the boundaries, the units of the string of sarcomeres have only one neighbor each, which leads to the following boundary conditions:

In additional simulations, we account for modifications of the shape of the cellular volumes of the compartments due to length changes of the sarcomeres. The cell length (longitudinal direction) was reformulated as a function of the fractional changes in sarcomere lengthwith the sarcomere resting lengths and SLi the sarcomere length for each unit i ∈ [1, 2, …, 50]. By maintaining the volume of the cell constant, we obtain the cell radius:with Vcell = 660.12 μm3.

As a result of these changes, the areas connecting the various compartments differ dependent on the changes in sarcomere length, and thus, the time-dependent diffusion parameters for each unit i ∈ [1, 2, …, 50] needed to be adjusted to reflect the changes in calcium diffusion between the compartments:with k ∈ [SS, cyto, SR]; τdyad,sarcol and τSS,cyto the original time-dependent diffusion parameter for diffusion within one unit from the dyadic cleft (dyad) to the SS and SS to the dyad, respectively; τSS, τcyto, and τSR the original time-dependent diffusion parameter for diffusion between adjacent SS, cytosolic, and SR compartments, respectively; and the updated values for the same time-dependent diffusion parameters for each unit i ∈ [1, 2, …, 50].

Incorporation of mechanical heterogeneities

We introduce mechanical heterogeneities by allowing sarcomeres within a single myocyte to have different resting lengths. To induce sarcomere heterogeneity, we deviated from the average resting length of 1.89 μm (34) by 5 and 10%, which are in the physiological range measured in myofibrils by Rassier et al. (40). Therefore, the resting length of the 50 sarcomeres were chosen to be numbers within the intervals [1.796, 1.985 μm] and [1.701, 2.079 μm], respectively. Using the rand function of MATLAB, which returns a normalized value between 0 and 1 drawn from a uniform distribution, we changed the range of the distribution to a new distribution within the desired interval [a, b] by multiplying each value by the width of the new range, (b–a), and then shifting every value by a.

In the following, we refer to simulations of myocytes with the same resting length for all sarcomeres as homogeneous (model (H)), myocytes with the same resting length for all but two sarcomeres at unit 15 and 32 as homogeneous without 2 (model (H2)), and myocytes with the same resting length for all but 2 × 5 sarcomeres at units 13–17 and 32–36 as homogeneous without 10 (model (H10)). The specific sarcomeres with length varying from the homogeneous string all have the same resting length, which differs by ±5 or ±10% from the homogeneous sarcomeres. Cases in which all sarcomere resting lengths are allowed to vary by 5% as 5% are referred to as inhomogeneous (model (IH5)), and myocytes with sarcomere resting lengths deviated by 10% as 10% inhomogeneous (model (IH10)). For models (H), (H2), and (H10), seven different models were implemented with varying sarcomere resting lengths (95, 97.5, 100, 102.5, 105, 107.5, and 110% of its original value), and 11 random variations were implemented for models (IH5) and (IH10). (The distribution of the sarcomere resting lengths for model (IH5) and (IH10) can be found in the Supporting Material.)

All units of a myocyte were simultaneously electrically stimulated with 5 ms pulses of 9.5 mV, which were added as a stimulation current to the transmembrane voltage at a cycle length of 600 ms for 160 beats to reach the limit cycle. The limit cycle was ascertained by comparing the upstroke, plateau, and repolarization phase of the last two APs by calculating the MSE over time between both curves until MSE < 0.01. For simulations at different extracellular calcium concentrations, the calcium concentration was changed first before running the models to the limit cycle. (To assess the limit cycle values, please refer to the Supporting Material.)

All simulations were run in MATLAB R2018a using ode15 s to solve the ODE system. For the simulations to reach limit cycle, a variable time step was used and set to 0.01 ms for all the last beat and all other simulations. (The code can be found at https://bitbucket.org/VTimmermann/biophys2019.git.)

Stretch protocols

In the 1D single-cell model, mechanical heterogeneities were investigated by subjecting the model to a wide range of stretch protocols varying in magnitude and stretch timing relative to cell activation. Myocyte stretch was modeled as an increase in sarcomere length by acute square waveform stretch experiments (22). All stretches (10, 12.5, and 15%) were imposed 10 ms after the last electrical stimulus and released at different timings, which varied in 10 ms steps in the interval of 300–2000 ms after the last stimulus.

Data analysis

The vulnerability window for myofilament-triggered waves and the wave velocities are presented as mean ± standard deviation calculated in MATLAB R2018a using the functions mean and std. The calcium wave velocities and suprathreshold DADs were analyzed using the one-way analysis of variance (anova1) with pairwise comparison of the group means (multicompare) preconditioned by the Bonferroni method in MATLAB R2018a. A data point was assumed to be an outlier if it did not lie within three standard deviations of the mean (three-sigma rule). Finally, we refer to a variable v as normalized when , for which the minimal and maximal values are calculated over all myocytes for a particular model.

Results

Strong electromechanical coupling promotes calcium waves and spontaneous activity

Experimental data suggests that under normal conditions, neither spontaneous (5, 6, 7) nor evoked (8, 9, 10, 11, 12), local calcium releases from the SR propagate (3). Hence, we increased extracellular calcium from control (1.8 mM) by 0.5 mM steps to simulate calcium overload until myofilament-triggered calcium waves occurred at Cao = 4 mM in model (IH5) at stretches of 15%. During calcium overload at Cao = 4 mM, 10% stretch did not induce myofilament-triggered waves in either of the models. In both models with sarcomeric dyssynchronicity, DADs occurred at 15% stretch, whereas only model (IH10) showed myofilament-triggered waves at 12.5% stretch. In model (H), 12.5 and 15% square waveform stretch experiments did not affect cytosolic calcium concentration for sarcomere resting lengths below 105% of its original value, whereas higher sarcomere resting lengths (≥1.9845 μm) triggered additional beats. Owing to the uniformity of the mechanics, either no calcium was spontaneously released from the SR or calcium was simultaneously released from the SR in all units of the myocyte, resulting in an additional beat.

In the heterogeneous models, spontaneous calcium release from the SR in one or multiple units can cause the propagation of a calcium wave, resulting in a calcium-mediated DAD, but only at specific stretch-release timings. In this section, we excluded the data from myocyte 4 of model (IH5) because for some stretch-release timings, an additional beat occurred. In model (IH10), the vulnerable window (see Fig. 4 A), in which stretch-release resulted in calcium-mediated DADs, was between 519 ± 29 and 777 ± 41 ms for 12.5% stretch and between 439 ± 20 and 900 ± 21 ms for 15% stretch. Although no myofilament-triggered waves were observed at 12.5% stretch for model (IH5), at 15% stretch, the vulnerable window was between 510 ± 35 and 815 ± 206 ms. In addition, the greatest vulnerable window (between 420 and 910 ms) occurred for model (H) when sarcomere resting length was increased by 5%. Thus, as shown in Fig. 4, an increase in stretch (green squares and red dots) or an increase in sarcomere heterogeneity (blue triangles and red dots) increased the susceptibility for myofilament-triggered calcium waves significantly.

Figure 4

Statistical representation of the lower (left) and upper (middle) bound for the vulnerability window in which myofilament-triggered calcium waves occur at Cao = 4 mM, as well as the timing of the fastest myofilament-triggered calcium wave (right) for model (IH5) (15% stretch, blue triangles), model (IH10) (15% stretch, red dots), and model (IH10) (12.5% stretch, green squares). To see this figure in color, go online.

It is worth noting that not only do sarcomere heterogeneity and stretch magnitude seem to modify calcium wave velocity but also that in the uncoupled model, no waves were observed for any model and any stretch protocol. These observations suggest that mechanoelectric feedback (MEF) may, under certain conditions, be able to modulate intracellular calcium wave dynamics. When the model was extended to incorporate changes in interunit distances, the susceptibility to myofilament-triggered calcium waves increased. However, halving the stretch magnitude resulted in similar results compared with the case without altered interunit distances (data not shown).

In addition to the vulnerable window of myofilament-triggered calcium waves, we also investigated the effect of the stretch-release timing on the calcium wave velocity. The velocities of the calcium waves in Fig. 4 B are given in , where cell activation duration was calculated as the delay between the first and last activated unit. The average calcium wave velocity of (IH10) at 12.5% was 240 ± 117 , compared to 431 ± 226 for model (IH10) at 15% and 441 ± 208 for model (IH5) at 15%. In contrast, the calcium waves at the lower and upper bound propagated with 176 ± 50 and 173 ± 48 , respectively, for model (IH10) at 12.5%; 186 ± 53 and 170 ± 56 , respectively, for model (IH10) at 15%; and 193 ± 82 and 180 ± 82 , respectively, for model (IH5) at 15%. As illustrated in Fig. 5 A, the wave velocity increased from the lower bound of the vulnerability window to the peak velocity at 628 ± 11 ms for (IH10) (green squares) at 12.5% stretch, 650 ± 10 ms for (IH10) at 15% stretch (red dots), and 653 ± 20 ms for (IH5) at 15% stretch (blue triangles) before decreasing to the upper bound.

Figure 5

(A) Calcium wave velocity increased from the lower bound of the vulnerability window to the peak velocity and decreased to the lower bound of the vulnerability window. (B) The maximal CaT at the wave increased from the lower bound of the vulnerability window to the peak velocity and decreased to the lower bound of the vulnerability window. (C) The maximal SR calcium concentration at the wave increased from the lower bound of the vulnerability window to the peak velocity and decreased to the lower bound of the vulnerability window. (D) Calcium wave velocity depended on the magnitude of the CaT. (E) The magnitude of the calcium-mediated DAD depended on the magnitude of the CaT. (F) The magnitude of the calcium-mediated DAD depended on the calcium wave speed. To see this figure in color, go online.

Similar to the calcium wave velocity (Fig. 5 A), for all models, the maximal CaT at the wave (Fig. 5 B) and the maximal SR calcium at the stretch-release (Fig. 5 C) increased from the lower bound of the vulnerability window to a peak at around 640 ms and decreased to the lower bound of the vulnerability window. Furthermore, our data suggest a dependency of the calcium wave velocity and the magnitude of the calcium-mediated DAD on the magnitude of the CaT at the wave (Fig. 5, D and E, respectively). Finally, the amplitude of the calcium-mediated DAD seems to be dependent on the calcium wave velocity. In summary, stretch-release timing affected the calcium handling of the mycyte resulting in calcium waves, whereas the magnitude and velocity of the calcium wave modified the amplitude of the occurring DAD.

For myocytes 1 and 11 of model (IH10) at 15% stretch, Fig. 6, A–D, respectively, show that the time course and magnitude of the myofilament-triggered calcium waves were dependent on the stretch-release timing. Fig. 6, A and C illustrate calcium waves at the lower bound of the vulnerability window, and Fig. 6, B and D show the corresponding calcium waves with the fastest velocity for the same cell. The difference in cell activation duration between the lower bound and the fastest wave for myocyte 1 was ∼25 ms. For this case, the location and number of origins for both the earliest and fastest wave remained the same, but the wave velocity varied by 8%. For myocyte 11, the cell activation duration decreased by 345 ms (from left to right), resulting in wave velocities of 225 and 999 , respectively. This change in wave velocity can be explained by the increase in wave origins and the recruitment of calcium sparks while the wave is propagating.

Figure 6

(A) Myofilament-triggered calcium wave of myocyte of model (IH10) at 15% stretch at the lower bound of the vulnerability window. (B) The fastest myofilament-triggered calcium wave of myocyte 3 of model (IH10) at 15% stretch with 8% increase in wave velocity is shown. (C) The myofilament-triggered calcium wave of myocyte 11 of model (IH10) is shown at 15% stretch at the lower bound of the vulnerability window. (D) The fastest myofilament-triggered calcium wave of myocyte 11 of model (IH10) is shown at 15% stretch with a 340% increase in calcium wave velocity due to multiple initiation sites and a higher CaT. To see this figure in color, go online.

Owing to the increase of SR calcium content as calcium was resequestered by SR calcium-ATPase (SERCA) in the interval between the lower bound and the fastest wave, the fastest calcium wave showed a higher CaT compared with all other calcium waves. Additionally, increased calcium content in the SR in all units may result in multiple origins of the myofilament-triggered calcium wave because the probability increases that the threshold for spontaneous calcium release from the SR is reached. At later stretch-release timings, the SR is already depleted of calcium.

Calcium concentration modifies mechanically triggered electrical instabilities

To investigate the role of extracellular calcium concentration, we assessed the vulnerable window for mechanically triggered electrical activity at Cao = 4.2 mM. For all simulations of model (IH10) and 9 of 11 simulations of model (IH5), low-magnitude DADs occurred after pacing was stopped (dashed black line in Fig. 7 A). Applying an acute square waveform stretch of 1% 10 ms after the last electrical stimulation and stretch-release 200 ms before the initiation of the spontaneous calcium wave, in four of the nine samples of model (IH5), no DAD occurred. For model (IH10), 1 of the 11 samples was excluded because a beat was triggered at 1% stretch.

Figure 7

(A) Simulations of myocyte 1 of model (IH5) at Cao = 4.2 mM for no stretch (dashed black lines), 1% stretch (red lines), and 2% stretch (blue lines) for membrane potential (top left), averaged whole-cell cytosolic calcium concentration (bottom left), averaged whole-cell SR calcium concentration (top middle), averaged whole-cell calcium buffering on TnC (bottom middle), averaged whole-cell normalized force (top right) (here, normalized force refers to the maximal normalized force given full activation as described in (34)), and the stretch protocol (bottom right). (B) Calcium waves for no stretch (left), 1% stretch (middle), and 2% stretch (right) are shown. The wave velocity increases by ∼280% from no stretch to 1% stretch and again by ∼180% from 1% stretch to 2% stretch. (C) Simulations of myocyte 1 of model (IH5) at Cao = 4 mM are shown for inhibited SERCA (blue lines), stimulated SERCA (red lines), and control (dashed black lines) for membrane potential (top left), averaged whole-cell cytosolic calcium concentration (bottom left), averaged whole-cell SR calcium concentration (top middle), averaged whole-cell calcium buffering on TnC (bottom middle), averaged whole-cell normalized force (top right) as described in (34), and the stretch protocol (bottom right). To see this figure in color, go online.

For the remaining five samples of model (IH5), we observed a difference in calcium wave magnitude by 33 ± 3.5%, in time to the peak of the calcium wave by 7 ± 3%, in force magnitude by 56 ± 3%, and in time to the peak of the force by 7 ± 6%. For the one sample of model (IH10), the time to the peak of the calcium wave varied by 2 ± 4%, the magnitude of the calcium wave by 36 ± 21%, the time to the peak of the force by 1.5 ± 3%, and the magnitude of the force by 31 ± 2%. The magnitude of the occurring DAD differed by 2.5 ± 3 and 5.4 ± 1.5% for models (IH5) and (IH10), respectively. The time to peak varied by 8 ± 3.5% for model (IH10) and decreased by 0.5 ± 1.5% for model (IH5).

When the stretch was increased from 1 to 2% (red and blue lines in Fig. 7 A, respectively), a beat was triggered, resulting from the increased magnitude and velocity of the myofilament-triggered wave. Fig. 7 B shows, from left to right, the increase in calcium wave velocity due to more synchronous SR calcium release, recruitment of sparks, and increased CaT magnitude. The increased CaT magnitude and wave velocity resulted in an increase of the calcium-mediated DAD amplitude. Consequently, myofilament-triggered calcium waves promoted a more synchronous calcium release, leading to faster wave propagation and potentially to suprathreshold DADs. Furthermore, stretch magnitude modified the velocity of the myofilament-triggered calcium waves (Fig. 7 B, 1% stretch (middle) and 2% stretch (right)), and myofilament-triggered waves propagated faster than spontaneous calcium waves because of the increase of SR calcium concentration at stretch-release.

It is worth noting that although for the calcium waves in Fig. 7, A and C, the SR and TnC calcium concentrations did not recover to control concentration, the concentration recovered when an additional beat occurred. Because of the low amplitude of the calcium-mediated depolarization resulting from the calcium wave, the membrane potential remained below the threshold of activating the voltage-sensitive L-type channels. Thus, less calcium entered the cell compared to stimulated beats. Nevertheless, after the spontaneous SR calcium release, calcium left the cell through the sarcolemmal exchange proteins, reducing the intracellular calcium concentration. Therefore, less calcium remained in the cell and was reuptaken into the SR.

To investigate whether calcium SR load could affect the initiation of myofilament-triggered calcium waves, we inhibited and stimulated SERCA by 25% (blue and red lines in Fig. 7 C, respectively) compared with control (dashed black lines). When SERCA was inhibited, spontaneous calcium release from the SR was prevented because of the reduced calcium content and no activation of neighboring release sites. On the other hand, consistent with experimental data (41), stimulation of SERCA increased calcium reuptake and stimulated neighboring calcium release sites by calcium diffusion into adjacent units with lower SR calcium concentration.

Additionally, we discovered that an increase of cytosolic calcium concentration at the end of a stimulated CaT while SR load was not fully recovered accelerated SR reuptake, resulting in a higher SR calcium concentration compared to control (see Fig. 7, A and C). During stretch, calcium buffering by TnC and SR calcium release was slightly increased, whereas the decline of the CaT was accelerated. After the twitch before stretch was released, the calcium concentration buffered on TnC was still elevated and became available when stretch was released. When the SR was not fully recovered, the resulting increase in cytosolic calcium load caused SR reuptake by SERCA, resulting in a further increase in SR calcium content. If the threshold for SR calcium to trigger activation of the RyR clusters was reached, opening of the RyR clusters was triggered, and calcium was released into the dyadic space. From the dyadic space, calcium first diffused into the SS and cytosol before propagating into the adjacent units of the myocyte.

To examine the dependence of the origins of the calcium waves on sarcomere heterogeneity, we ran square waveform stretch experiments for models (H2) and (H10). Although the behavior of model (H2) was similar to model (H), for model (H10), calcium waves occurred when the majority of the sarcomere resting lengths were increased by 10% of the original value. For simulations at higher resting lengths, an additional beat was triggered. For the groups with a lower sarcomere length (5% of the original value) compared to the majority, the waves originated from the middle of these two groups. The origins were located at the units with the maximal sarcomere length change (red lines in Fig. 8 A) and the highest SR calcium concentration (blue lines in Fig. 8 A). On the other hand, for simulations in which the two groups had a longer sarcomere resting length, the origins of the calcium wave were located outside of the groups. Yet the calcium waves originated from the units with the maximal sarcomere length change (red lines in Fig. 8 C) and the highest SR calcium concentration (blue lines in Fig. 8 C).

Figure 8

(A) Maximal SR calcium concentration (blue lines) and maximal sarcomere length changes (red lines) for each unit of the myocyte for simulations of model (H10) with two groups of sarcomeres with a shorter resting length (1.7955 μm) compared to the resting length of all other sarcomeres in the myocyte (1.89 μm). (B) The calcium wave of model (H10) with two groups of sarcomeres with a shorter resting length (1.7955 μm) is compared to the resting length of all other sarcomeres in the myocyte (1.89 μm). (C) Maximal SR calcium concentration (green lines) and maximal sarcomere length changes (red lines) for each unit of the myocyte for simulations of model (H10) with two groups of sarcomeres with a longer resting length (1.9845 μm) are compared to the resting length of all other sarcomeres in the myocyte (1.89 μm). (D) The calcium wave of model (H10) with two groups of sarcomeres with a longer resting length (1.9845 μm) is compared to the resting length of all other sarcomeres in the myocyte (1.89 μm). To see this figure in color, go online.

Our results suggest that SR calcium concentration caused by calcium reuptake by SERCA during stretch-release is the main determinant of accelerated velocities of myofilament-triggered calcium waves (Fig. 9). The calcium wave velocity changes observed in our simulations can all be explained by varying SR calcium concentrations caused by differences in calcium buffering to TnC, which were in turn caused by varying sarcomere lengths. This result is consistent with our observation that at increased resting sarcomere lengths, spontaneous beats occurred as the result of synchronous SR calcium release in all units of the cardiomyocyte. Compared to heterogeneous simulations, calcium was released at once from all units of the myocyte when the threshold for calcium SR content for spontaneous release was reached. Owing to the uniformity of these events, not only were additional beats triggered, but the vulnerable window also increased.

Figure 9

Schematic illustration of myofilament-triggered calcium waves. The red arrow indicates calcium release from the myofilaments due to stretch-release. The black arrows show the pathway for calcium reuptake to the SR via SERCA. If the resulting increase in calcium content in the SR reaches the threshold for spontaneous calcium release, the RyR clusters are activated and release calcium into the dyadic space. The orange arrow shows the diffusion pathway of calcium from the dyadic space into the cytosol of the same unit and then into the adjacent unit. In the neighboring unit, calcium release from the SR may again be triggered by an increase in cytosolic calcium concentration, which results in calcium reuptake to the SR. RyR, ryanodine receptor. To see this figure in color, go online.

Discussion

The presented simulation study demonstrates that dynamic mechanical coupling in spatially explicit myocyte models with mechanical perturbations may influence calcium wave velocity and trigger supratheshold DADs that may be arrhythmogenic. In particular, stretch and release of mechanically heterogeneous myocytes may affect vulnerability to arrhythmogenic calcium waves. Associated with stretch and release of sarcomeres in dynamically coupled EM, strain-rate dependent changes in calcium sensitivity of the myofilaments can contribute to calcium release from TnC to trigger calcium release from the SR, resulting in an intracellular calcium wave through the combined effect of time-dependent calcium diffusion and CICR. Further, intracellular heterogeneities in mechanics may promote additional instability through increased wave velocity, but only at particular stretch timings and during calcium overload.

Contractile propagation process and wave velocity

Experiments have shown that calcium waves originate close to border zones of mechanically heterogeneous tissues and can propagate through the tissue to cause premature beats (1,23). These findings suggest that MEF, and in particular mechanics-induced calcium waves, is important to fully understand the mechanisms leading to arrhythmia resulting from mechanical heterogeneities (26,30,42, 43, 44). This study examines myofilament-triggered calcium waves in single cells and how these may be initiated by cell stretch and subsequent release, applied after electrical stimulus of the cell. ter Keurs, Boyden, and others (25,45,46) have pointed out that abnormal stretch patterns could be expected in regions of myocardial mechanical heterogeneity—for instance, in the border zone of ischemic regions—and may provide a trigger for calcium waves and DADs that originate from these regions. Our results also suggest that the velocity of calcium wave propagation affects the DADs’ amplitude, consistent with experimental observations (26,30,42,47).

Spontaneous calcium release events are associated with a calcium-loaded SR (3,5,13,18, 19, 20), even though they were first observed as spontaneous contractile waves (16,17,21,48). Experiments by Cheng et al. (3) indicate that calcium sparks occur at the site of wave initiation, similar to our finding of spontaneous calcium waves at Cao = 4.2 mM. These slowly propagating calcium waves are characterized by a discrete calcium release from the SR forming the wave front, recruiting other sparks along the way by progressing in a saltatory manner, as has been shown before by (3). Fast-propagating calcium waves occurred when the cell was exposed to rapid changes in length, consistent with previous studies (22,23,26,29,30). In addition to the known correlation between extracellular calcium concentration and wave velocity (3,5,13,18, 19, 20), our study suggests that mechanical heterogeneity at the sarcomeric level and sarcomere length distribution can modulate wave speed.

Experiments in inhomogeneous cardiac muscle discovered arrhythmogenic contractile waves, suggesting that the arrangement of weaker damaged regions in series with normal strong muscle units is responsible for spontaneous contractile waves (1). The underlying mechanism has been hypothesized to be caused by calcium release: first by the myofilaments because of stretch and release, and then by the SR in the cells near the damaged region. Hence, tissue-level heterogeneities in mechanical loading may trigger spontaneous contractile waves in stretched myocytes, which may lead to propagating CaTs that can be conducted into neighboring myocytes (22,27,29).

Using single-cell preparations, experimental studies have characterized calcium oscillations and contractile waves in cardiomyocytes (16,49,50). Computational efforts to describe these spatiotemporal phenomena in a mechanistic manner have included simulations of EP in elementary and more detailed cellular geometries (51, 52, 53). But most of these mathematical models did not incorporate mechanics or MEF. Integrating experimental data on mechanical and electrophysiological processes with appropriate constitutive models provides a framework to analyze MEF in the propagation of calcium waves. Recent efforts proposed schemes for contractile cardiac cells in continuum mechanics frameworks (54,55). Others described calcium sparks and waves in detailed ventricular cell models of spatially explicit distribution of RyR clusters (56,57). Finally, to include components of the electrophysiological response to cardiac cell contractility, single strain-energy functions have been integrated to computational models (58). Focusing only on the longitudinal propagation of calcium waves in one dimension results in a lack of information about cellular anisotropy and radial wave propagation (35). However, two-dimensional and three-dimensional models did not include MEF or mechanically triggered calcium waves.

Our results for single-myocyte simulations demonstrated that in the presence of subcellular mechanical heterogeneities, abnormal stretch and release may trigger spontaneous calcium waves. These waves may in turn trigger contractile waves and suprathreshold DADs because of their fast propagation velocity. Abnormal patterns of cell stretch may be present in regions of mechanical heterogeneity, for instance, during acute myocardial ischemia. Therefore, as previous investigators have hypothesized, MEF may contribute to coupling regional mechanical dysfunction with the initiation of premature beats and arrhythmias in the adjacent myocardium (25,45,46).

Limitations

Although this study provides a range of, to our knowledge, new insights regarding MEF arrhythmogenesis, our approach has some limitations. A significant limitation was our choice of the contractile myofilament model (34). Although this mechanical model is computationally fast, it relies on a simplified representation of calcium binding to the myofilaments and is primarily built using steady-state data. More recent mechanical myofilament models include 1) an explicit spatial modeling of nearest-neighbor interactions within an ODE-based model framework (59), which represents a better structure for 1D model extensions; 2) calcium-independent activation of the thin filament in a Monte Carlo model with nearest-neighbor interactions (60), which allows for stronger MEF; 3) thin and thick filament compliance, including XB-XB cooperativity and three-dimensional spatial elements with lattice spacing, which provides a configuration for strain development resulting from changes in lattice spacing (61); and 4) an explicit formulation for thin and thick filament kinetics incorporating a hypothesized super-relaxed state (62).

The EP model, on the other hand, may reduce the ability of mechanically induced changes in CaTs to affect RyR opening because of its formulation of the SS. Owing to the current formulation of the diffusion between the SS and cytosol, as well as the calcium concentration in those two compartments, calcium diffusion from the cytosol into the SS is unlikely. Furthermore, the chosen EP model produces an increase in AP duration with increasing extracellular calcium. This has been noted as a current limitation to numerous modern models (63), and because our simulations were run at elevated calcium, this could have an important effect on the vulnerability window of calcium-driven DADs. Additionally, in the Shannon et al. model (32), calcium diffusion from the cytosol into the SS is relatively slow because the original model is constructed to represent the higher calcium concentrations near the sarcolemma. But the Shannon et al. model (32) was not built to reproduce MEF mechanisms, and thus, diffusion from the cytosol into the SS was not the main focus of the original model. Therefore, the initiation of calcium waves more likely occurs as a result of spontaneous calcium release from the SR due to calcium reuptake via SERCA than CICR through the activation of the RyRs. This limitation is also responsible for the delay between the moment of stretch-release and the occurrence of the DAD.

In addition to myofilament regulation of calcium, SACs are known to contribute to mechanics-induced arrhythmia by depolarizing the diastolic membrane potential in myocytes (and whole hearts) (2,64, 65, 66). However, proarrhythmic effects of stretch have also been reported in experiments with chemically blocked SACs (22,25). Also, in addition to SACs, another known mechanism that contributes to MEF are alterations in the open probability of intracellular calcium release mechanisms, resulting from signaling processes collectively known as X-ROS (67,68). In this study, X-ROS has not explicitly been simulated because we used whole muscle data for the parameterization of myofilament-triggered calcium release. It is likely that X-ROS-mediated responses are included within those dynamics; however, it is not possible for us to separately analyze these dynamics, and it is likely that X-ROS is more important than myofilament calcium sensitivity in some contexts, particularly diastolic stretch (67,68).

Another limitation of this study is that we did not test more physiological stretch patterns to investigate how experimentally observed stretches in disease effect myofilament-triggered calcium release and AP duration. It is possible that for less rapid stretch patterns, spontaneous beats could occur, or a spontaneous beat could be triggered when stretch is released during a DAD in a shorter vulnerability window. Additionally, we did not impose different pacing rates, even though experimental data (16) suggest that with increased pacing rate, the probability of spontaneous mechanical oscillations occurring simultaneously at multiple foci increases. These contractile oscillations could partially synchronize, resulting in contractile waves with increased amplitude. Moreover, our conclusions were based on a deterministic model that describes the average behavior of a process. But low intracellular calcium concentrations and the small volumes of intracellular subdomains mean that the number of calcium ions in the cleft is low and intrinsically stochastic such that the calcium-handling processes can vary notably from the average in single simulations. Mathematical simulations that account for this stochasticity have found that for small compartments such as the dyadic cleft, the opening probability of receptors such as the RyRs depend on whether calcium ions are present and not if those calcium ions are in close proximity. Hake and Lines (69) showed that these notable variations between the stochastic and deterministic calcium diffusion processes vanish on larger timescales such as the ones presented in this study. Therefore, stochastic approaches would be useful for investigating the effects of mechanical heterogeneity, though this would probably also require considering more strain nonuniformity on a more finely resolved submicron scale than we have considered here. Nevertheless, the lack of a stochastic model of myofilament-triggered calcium waves remains a significant limitation, especially in the comparisons we have made between homogeneous dynamics in which we cannot reproduce propagating waves that do not activate the entire cell and heterogeneous systems that lack spontaneous RyR opening that is not triggered by increased SR calcium load. This should be investigated in future studies to ensure that our main conclusions hold in the presence of stochastic RyR gating.

Finally, the origin of calcium waves in experiments of injured sarcomeres in series with normal muscle tissue occur in proximity to the injury border zone. Because this study only investigates single cells, the complete mechanism of these mechanically driven tissue-scale calcium waves could not be captured and needs further investigation.

Conclusions

In this study, we investigated the role of subcellular mechanics underlying changes of mechanically triggered calcium waves and how mechanical heterogeneities may modify wave velocity compared to spark-initiated waves. The goal was to investigate whether the changes in myofilament calcium sensitivity contribute to proarrhythmic effects in calcium overload conditions. Using strongly coupled EM models of rabbit ventricular myocytes, we showed that the myofilament-triggered waves can propagate with higher velocities than spark-initiated waves, but only with heterogeneous mechanics and at specific stretch-release timings.

The simulation results indicate that in conditions of calcium overload, an increase in stretch or an increase in sarcomere heterogeneity increases the electrical susceptibility for myofilament-triggered waves significantly depending on the specific stretch-release timing. Furthermore, the specific sarcomere heterogeneity can modulate the calcium wave susceptibility in the different units of the myocyte and affects aftercontraction wave propagation velocity. These results indicate that strongly coupled EM may play a role in modulating calcium dynamics in time and space during physiologically plausible dynamic stretch.

Author Contributions

V.T. carried out all simulations, analyzed the data, and wrote the article. A.G.E. commented and edited the manuscript. S.T.W. discussed the methods and helped interpreting the results. J.S. discussed the methods, helped interpreting the results, and made comments to the manuscript. A.D.M. designed the research and commented and edited the manuscript.