It is not understood why, after onset of left bundle-branch block (LBBB), acute worsening of cardiac function is followed by a further gradual deterioration of function, whereas most adverse cardiac events lead to compensatory adaptations. We investigated whether mechano-electrical coupling (MEC) can explain long-term remodeling with LBBB and cardiac resynchronization therapy (CRT). To this purpose, we used an integrative modeling approach relating local ventricular electrophysiology, calcium handling, and excitation-contraction coupling to global cardiovascular mechanics and hemodynamics. Each ventricular wall was composed of multiple mechanically and electrically coupled myocardial segments. MEC was incorporated by allowing adaptation of L-type Ca2+ current aiming at minimal dispersion of local external work, an approach that we previously applied to replicate T-wave memory in a synchronous heart after a period of asynchronous activation. LBBB instantaneously decreased left-ventricular stroke work and increased end-diastolic volume. During sustained LBBB, MEC reduced intraventricular dispersion of mechanical workload and repolarization. However, MEC-induced reduction in contractility in late-activated regions was larger than the contractility increase in early-activated regions, resulting in further decrease of stroke work and increase of end-diastolic volume. Upon the start of CRT, stroke work increased despite a wider dispersion of mechanical workload. During sustained CRT, MEC-induced reduction in dispersion of workload and repolarization coincided with a further reduction in end-diastolic volume. In conclusion, MEC may represent a useful framework for better understanding the long-term changes in cardiac electrophysiology and contraction following LBBB as well as CRT.
- cardiac resynchronization therapy
- excitation-contraction coupling
- left bundle-branch block
- mechano-electrical coupling
cardiac resynchronization therapy (CRT) has emerged as an important therapy to improve pump function in patients with chronic heart failure and conduction disturbances, such as left bundle-branch block (LBBB) (10). This therapy aims to restore synchrony of ventricular electrical excitation, generally by pacing the left ventricle (LV) at a single location (LV pacing) or by pacing both ventricles (simultaneous or sequential biventricular pacing). A common observation is that the improvement in pump function when starting CRT, for example, quantified as increase in LV ejection fraction, further continues over time (8). A similar, but reverse pattern is seen after stopping CRT in these patients (61).
In the canine model of LBBB (58), LBBB resulted in an acute reduction of LV pump function, followed by a further deterioration over time (59), whereas CRT caused an acute improvement, followed by ongoing improvement in the long run (58). Such behavior seems different from other cardiovascular abnormalities, like hypertension or infarction, where cardiac adaptation, at least initially, leads to compensation of the functional impairment (36, 39).
Various experimental studies have shown extensive and differential changes in gene and protein expression as well as posttranslational changes in dyssynchronous hearts (2, 9, 35, 51), resulting in abnormal calcium handling and repolarization. These findings do not, however, clarify directly why cardiac pump function would worsen in chronic LBBB.
In the present study, we investigated the relation between local electrophysiological and global functional alterations induced by LBBB and CRT using an integrative modeling approach. To that end, we extended our electromechanical model of the human left heart and systemic circulation (20) with the right heart and pulmonary circulation. The present model also includes mechanical interaction of the LV and right ventricle (RV) through the interventricular septum as well as hemodynamic interaction between the two ventricles through the systemic and pulmonary circulations (33, 34). In a previous study, we showed that mechano-electrical coupling (MEC), i.e., remodeling of ionic currents and calcium handling triggered by changes in local mechanical workload, explains optimal tuning of cardiac pump function as well as the concordant T-wave during normal sinus rhythm and T-wave memory after temporary abnormal electrical activation (20). In particular, it was assumed that specific conductance of the L-type calcium current (GCaL) increases when mechanical load is too low and decreases when workload is too high, aiming at minimal dispersion of local mechanical external work. Our choice for changing GCaL in that study and in the present study was based on the observation that expression of L-type calcium channels and the calcium transient show regional differences in dyssynchronous heart failure (DHF) attributable to LBBB and rapid pacing in dog hearts, complying with local regulation (2), and that calcium channel blockade attenuates cardiac memory after temporary pacing (40). The aim of the present study was to investigate whether MEC by means of local regulation of L-type calcium current conductance could explain the deterioration of function during LBBB as well as the long-term improvement during CRT.
MATERIALS AND METHODS
The model used in the present study is based on our CircAdapt model of human cardiovascular mechanics and hemodynamics (6, 7, 33), incorporating mechanical interaction of the LV and RV through the interventricular septum and hemodynamic interaction between both ventricles through the systemic and pulmonary circulation. Our model has a modular setup, describing all elements of both circulations including valves, arteries, microcirculation, veins, and atria. Global ventricular pump mechanics (pressure-volume relation) was related to myofiber mechanics (myofiber stress-strain relation) using the principle of conservation of energy (4, 5). To simulate fast processes such as ionic channel dynamics, calcium handling, and excitation-contraction coupling combined with the much slower process of cardiac electromechanical remodeling, we have adopted a geometrically simplified model of the heart. The model includes two atrial walls and three ventricular walls: LV free wall (LVfw), interventricular septum, and RV free wall (RVfw) (33). Electromechanical behavior of a single ventricular wall is modeled by a multitude of mechanically and electrically coupled myocardial segments (Fig. 1). The number of segments for each ventricular wall was proportional to its wall volume (Vwall): 120 segments for the RVfw (Vwall,RVfw = 42 ml), 150 segments for the septum (Vwall,septum = 50 ml), and 300 segments for the LVfw (Vwall,LVfw = 101 ml).
In the present study, the heuristic model of cardiac mechanics in the original CircAdapt model (33) was replaced by a cascade of physiological models (Fig. 1). Ionic currents and calcium handling were modeled by the 2006 model of the human ventricular action potential from ten Tusscher et al. (54, 55). Mechanical behavior of a segment was modeled by the classical three-element rheological scheme (21, 49). Excitation-contraction coupling was modeled by model 5 of Rice et al. (46) and depended on the intracellular concentration of free calcium [from the model of ten Tusscher et al. (55)], sarcomere length, and velocity of sarcomere shortening. These models were integrated to describe electromechanical behavior of multiple segments that were electrically and mechanically coupled in series as described previously (20, 25, 28, 29, 30). Intracellular and extracellular electrical coupling between the segments was modeled according to the bidomain equations (19, 27).
Cardiac mechanics and hemodynamics.
The geometry of the biventricular model was described by wall volume (Vwall) and reference surface area of the midwall (Am,ref) for LVfw, septum, and RVfw, where Am,ref is defined as the quotient of Vwall and the wall thickness at mid-ejection. Parameters describing ventricular geometry were based on human data. Wall volumes were obtained from Doherty et al. (15). Am,ref was found assuming wall thickness at mid-ejection of 1.3 cm for LVfw (52), 1.1 cm for septum (52), and 0.4 cm for RVfw (15) (Table 1).
Adaptation of electromechanical properties (in the present study, specific conductance of L-type Ca2+ current) requires proper boundary conditions. In our model, these boundary conditions were imposed by pressure control and volume control. With pressure control, the systemic peripheral resistance was adjusted to maintain a predefined mean arterial pressure (MAP). Also the amount of circulating blood was regulated to obtain the required cardiac output (CO). Adaptation of circulating blood volume represented the combination of short-term blood flow regulation (to control venous return) and long-term volume control by the kidneys. Using these regulatory mechanisms, CO and MAP were maintained throughout the series of simulations. In the present study, CO was 5.1 l/min (85 ml/s), MAP was 91.5 mmHg (12.2 kPa), and cardiac cycle time (tcycle) was 0.85 s (heart rate 71 beats/min). These settings were the same as in Hermeling et al. (20) and in Lumens et al. (33).
Ventricular electrophysiology was modeled using our previously published bidomain model (26, 27). As described above, electromechanical behavior of a ventricular wall was modeled by a serial arrangement of mechanically and electrically coupled myocardial segments. All segments within a wall could be either activated simultaneously or as a consequence of simulated impulse propagation. In the second case, activation time of a segment depended on the distance to the site of electrical stimulation and the bidomain parameters (Table 1).
The electrophysiological state of each segment was defined by the intracellular potential (Vint), the extracellular potential (Vext), and the state of the cell membrane, which was expressed in gating variables and ion concentrations. The membrane potential was defined by Vmem = Vint − Vext. Exchange of current between the intracellular and extracellular domains occurred as transmembrane current (Itrans), which consisted of ionic currents (Iion) and capacitive currents: (1) where χ is the ratio of membrane surface to tissue volume and Cmem is membrane capacitance per unit of membrane surface (Table 1). Intracellular and extracellular currents between neighboring segments were related to intracellular and extracellular conductivities (gint and gext in Table 1).
To model ionic currents and calcium handling in ventricular tissue, we applied the 2006 model of ten Tusscher et al. (54, 55). The total ionic current was given by (2) where INa is fast Na+ current, IK1 is inward rectifier K+ current, Ito is transient outward K+ current, IKr is rapid delayed rectifier K+ current, IKs is slow delayed rectifier K+ current, ICaL is L-type Ca2+ current, INaCa is Na+/Ca2+ exchanger (NCX) current, INaK is Na+/K+ pump current, IpCa is Ca2+ pump current, IpK is K+ pump current, and IbCa and IbNa are background Ca2+ and Na+ currents (54, 55). The main difference between the 2004 and 2006 model of ten Tusscher et al. (55, 56) lies in the formulation of calcium dynamics and ICaL. The 2006 model includes subspace calcium dynamics that controls ICaL and calcium-induced calcium release (CICR), modeled with a four-state Markov model for the ryanodine receptor (RyR). In addition, both fast and slow voltage-gated inactivation of ICaL are included (55).
For each ventricular wall (RVfw, septum, LVfw), electromechanical behavior of the midwall was assumed. The midwall was defined as the surface through the wall such that wall mass on the endocardial site was equal to wall mass on the epicardial site. The model of ten Tusscher et al. (55) distinguishes between endocardial, midwall, and epicardial cells, which differ in Ito and IKs. In the normal canine heart, there seems to be a gradual transmural change in electrophysiological properties from endocardium to epicardium (31), which may serve to obtain synchronized contraction across the ventricular wall (13). Because the midwall in our model is more located toward the epicardium, we have chosen to model cellular electrophysiology in between the midwall cell and the epicardial cell (75% midwall cell and 25% epicardial cell).
In the model of ten Tusscher et al. (54, 55), the epicardial, midwall, and endocardial cells differ in specific conductance of IKs an Ito (GKs and Gto, respectively). In the present study, Gto = 0.294 nS/pF is equal to the value used by ten Tusscher et al. (55) for epicardial and midwall cells, and GKs = 0.172 nS/pF is in between the value for epicardial cells (0.392 nS/pF) and midwall cells (0.098 nS/pF) (55). To investigate the sensitivity of our model with respect to variation in specific conductance of potassium currents, we also performed a series of simulations in which Gto and GKs were equal to values chosen for endocardial, midwall, and epicardial cells. Computational aspects of electrophysiology are discussed in appendix a.
Electromechanical behavior of a ventricular wall was modeled by a serial arrangement of mechanically and electrically coupled myofiber segments. As previously described (20, 28, 29, 30), mechanical behavior of a segment was modeled by the classical three-element rheological scheme (21, 49) composed of a serial arrangement of a contractile (CE) and a series elastic element (SE), placed in parallel to a parallel elastic element (PE). Forces generated by each of these elements are expressed as first Piola-Kirchhoff stress, denoted by TCE, TSE, and TPE, respectively. The length of the PE is the sarcomere length and is denoted by ls.
Contractile stress (TCE) generated by the CE depended on intracellular Ca2+ concentration ([Ca2+]i), internal sarcomere length (lsi, defined as the length of the CE), and velocity of sarcomere shortening v = −dlsi/dt. TCE was defined by (3) where fCE is a scaling factor, fv(v) is Hill's force-velocity relation, and Fnorm is normalized contractile force generated by the sarcomeres. Function fv(v) was defined by (4) where vmax is maximum velocity of sarcomere shortening and cv is a constant describing the shape of the hyperbolic relationship (22) (Table 1). Fnorm, related to [Ca2+]i and lsi, was modeled by model 5 from Rice et al. (46) with parameters obtained from Rice et al. (44).
First Piola-Kirchhoff stress present in the series elastic element (TSE) was determined by the average cross-bridge length, lSE = ls − lsi, and the number of active cross-bridges, which was assumed to be proportional to Fnorm: (5) where fSE is a scaling factor and kSE a material constant (Table 1).
Stress generated by the PE (TPE) was the summation of stress generated by titin (TPE,titin) and stress generated by collagen (TPE,collagen), which were defined by (6) (7) where fPE is a scaling factor, kPE,titin, and kPE,collagen are material constants describing the elasticity of titin and collagen, and ls0 and lPE0,collagen are reference lengths (Table 1). First Piola-Kirchhoff stress generated by the segment (Tsegment) was the sum of stress generated in the parallel elastic and in the contractile elements. In summary, it holds for the three-element model: (8) (9)
A computational scheme to solve the equations describing myofiber mechanics is derived in appendix b.
MEC was applied by increasing or decreasing specific conductance of the L-type Ca2+ current (GCaL) when external work (Wext) was below or above the target value, respectively. Adaptation of GCaL was done automatically during the simulations. After each cardiac cycle, Wext, defined as the area of the local stress-strain loop, was computed for each segment, and GCaL was adjusted. In case segmental Wext was below the target value, GCaL was augmented, allowing increase of L-type Ca2+ current (ICaL) for that segment. In case Wext was above the target value, GCaL was decreased. The target value for Wext was based on the assumption that the ventricles in total generate 1.5 J work per beat. With a total ventricular wall volume of 193 ml (Table 1), this leads to target external work (Wext,ref) of 7.8 kJ/m3.
In the present study, GCaL was allowed to deviate ±20% from the default value. This range was based on experimental results from Aiba et al. (2), who measured ICaL in control dogs and in dogs with DHF. They measured an increase in ICaL amplitude in the LV anterior wall (−4.0 ± 0.2 pA/pF in control and −4.6 ± 0.4 pA/pF in DHF) and a decrease in ICaL amplitude in the LV lateral wall (−3.8 ± 0.1 pA/pF in control and −3.2 ± 0.2 pA/pF in DHF).
Adaptation of GCaL was implemented as follows. Initially, GCaL,n = GCaL,ref for each segment n, with GCaL,ref the reference value (55). Each time a new cardiac cycle started, GCaL,n was adapted as follows: (10)
To ensure that GCaL,n was in steady state for each segment n, 400 cardiac cycles were simulated. In the steady state, each segment n either had reached the target value Wext,ref or had a GCaL,n that had reached its lower or upper limit (0.8·GCaL,ref and 1.2·GCaL,ref, respectively).
In Fig. 2, the effect of varying GCaL on action potential morphology, calcium transient, and Fnorm is shown for endocardial, midwall, and epicardial cells as modeled by ten Tusscher et al. (55) as well as for our settings of GKs and Gto. We set GCaL to 80%, 100%, and 120% of the default value, and lsi was set to 2.2 μm. Action potentials were generated by application of a stimulus current of −38 pA/pF during 1 ms with a stimulus interval of 1,000 ms. Data are shown for the fiftieth action potential.
To investigate the effect of remodeling of GCaL on local electromechanical function as well as on global LV pump function, a series of simulations was performed consecutively: 1) normal sinus rhythm, 2) acute LBBB (without MEC), 3) sustained LBBB (with MEC), 4) acute CRT (without MEC), and 5) sustained CRT (with MEC). During normal sinus rhythm and acute LBBB, all segments had the default value for GCaL (55). Then, GCaL was adapted by MEC during sustained LBBB. These new settings for GCaL were then retained during acute CRT because, in clinical practice, CRT is usually applied in patients with chronic LBBB. Finally, GCaL was adapted by MEC during sustained CRT.
Normal sinus rhythm was simulated by simultaneous activation of all segments at 0 ms, i.e., 155 ms after stimulation of the right atrium (Fig. 3A). Abnormal impulse propagation, as occurs in LBBB, was simulated by electrically stimulating the middle of the septum and scaling conductivity parameters so that it took 40 ms for the action potential to propagate to the junction with the RVfw and the LVfw and another 80 ms to activate the most remote part of the LVfw. As a consequence, total activation time increased to 120 ms. This degree of dyssynchrony is similar to that used in a previous simulation study by Lumens et al. (34), albeit with a small time delay between RVfw and septum. Moreover, the 120 ms in the model of transmurally averaged activation times matches with the QRS duration of ∼150 ms in human with complete LBBB (53). During LBBB, RVfw electrical activation was simulated as in normal sinus rhythm (Fig. 3B).
Biventricular pacing during CRT was modeled by simultaneous stimulation of the septum and the lateral site of the LVfw using the same conductivity parameters as during LBBB. To model prolonged duration of RVfw activation during biventricular pacing compared with sinus rhythm (57), impulse propagation was simulated by stimulating the lateral site of the RV free wall. Using the same conductivity parameters, total RVfw activation time was 32 ms (Fig. 3C).
As an example of the effect of MEC on cellular behavior, Fig. 4 shows membrane potential (Vmem), calcium transient ([Ca2+]i), and normalized contractile force (Fnorm) for the center segments of septum and LVfw during acute and sustained LBBB. During acute LBBB, action potentials and calcium transients of the septal and the LVfw segments were similar, but, because of the delayed electrical activation, the LVfw segment was more prestretched, leading to a higher Fnorm (Frank-Starling mechanism). During sustained LBBB, MEC reduced the intraventricular differences in external work by increasing GCaL for the septal segment and decreasing GCaL for the LVfw segment. This process led to an increased Ca2+ transient and Fnorm in the septal segment and a decreased Ca2+ transient and Fnorm in the LVfw segment. In addition, action potential duration (APD-60mV) increased by 20 ms in the septal segment and decreased by 27 ms in the LVfw segment, such that repolarization of the two segments occurred more synchronously in the sustained LBBB simulation.
Stress and strain.
Figure 5 shows natural strain and Cauchy stress–natural strain relations, averaged over the septum and over the LVfw. During normal sinus rhythm, shape and area of the septal and LVfw stress-strain loops were similar, i.e., nearly the same amount of Wext was generated. During acute LBBB, the stress-strain loop area became larger in the late-activated LVfw and negative in the early-activated septum (figure-of-eight loop). During sustained LBBB, MEC led to an increase in stress-strain loop area for the septum. However, Wext in the septum was still smaller than normal. Moreover, after the initial increase of the loop area of the LVfw with acute LBBB, a decrease was induced by MEC. Upon starting CRT, Wext was larger in the septum than in the LVfw, but this difference nearly vanished during sustained CRT.
Effect of MEC during LBBB and CRT.
Figure 6 displays the distribution of Wext, GCaL, APD-60mV, and time of repolarization (trepol) during LBBB (left) and CRT (right) across the LV myocardium. With acute LBBB, Wext decreased in all septal segments and increased in all LVfw segments to a degree that was related to activation time. During sustained LBBB, MEC increased Wext in the septum, but for most septal segments target value Wext,ref was not reached because the maximum allowed increase of GCaL (120% of baseline) was not sufficient to overcome the effect of early onset of contraction. On the other hand, MEC was able to reduce Wext to the target value for all LVfw segments. As a consequence, average Wext generated by all LV segments together reduced from 8.4 kJ/m3 to 7.2 kJ/m3. With onset of CRT, Wext increased in most septal segments, whereas it decreased in most LVfw segments, but average Wext remained 7.2 kJ/m3. During sustained CRT, MEC led to normal Wext for all segments of the septum and the LVfw, except for the early-activated segments in the LVfw. This led to an increase of average work to 7.3 kJ/m3.
MEC-induced changes in GCaL affected APD-60mV, causing a reduction of APD-60mV in the late-activated LVfw and an increased APD-60mV in the early-activated septum during sustained LBBB. These changes consequently led to a smaller dispersion of repolarization during sustained LBBB compared with acute LBBB. Similarly, APD-60mV increased near the pacing sites during sustained CRT, also reducing dispersion of repolarization.
Figure 7 shows tracings of aortic and LV pressure, LV volume, and aortic flow velocity, as calculated by the model. With acute LBBB, systolic LV and aortic pressure decreased from 129 to 126 mmHg and peak aortic flow velocity from 1.22 to 1.12 m/s, whereas LV end-diastolic volume (EDV) increased from 132 ml to 137 ml. With sustained LBBB, systolic pressure further decreased to 121 mmHg and peak aortic flow velocity to 0.88 m/s, whereas EDV further increased to 153 ml, indicating further deterioration of LV function. Acute CRT increased systolic pressure to 124 mmHg and peak aortic flow velocity to 1.04 m/s, whereas EDV remained at 152 ml. Sustained CRT had no further effect on systolic pressure and aortic flow velocity. However, EDV decreased from 152 ml to 148 ml, which indicates improvement in function.
Figure 8 summarizes some of the previously mentioned results, depicting that MEC reduced the dispersion in Wext and in trepol by ∼50% during sustained LBBB and during sustained CRT compared with the respective acute phases. With LBBB, MEC led to a further deterioration of LV function as indicated by a decrease in stroke work and an increase in EDV compared with acute LBBB. Following the acute hemodynamic improvement by CRT, MEC further improved LV function, as evidenced especially by the decrease in EDV compared with acute CRT.
Figure 9 illustrates the contribution of the three ventricular walls to the total amount of Wext generated by the ventricles. During normal sinus rhythm, the differences in Wext generated by the ventricular walls reflect the differences in wall volume (42, 50, and 101 ml for RVfw, septum, and LVfw, respectively). During acute LBBB, Wext increased for the LVfw, whereas it decreased to zero in the septum. Note that Wext in the RV wall decreased, whereas total ventricular work increased, indicating that contraction in the LVfw also supports RV function. The increase in Wext for the LVfw diminished during sustained LBBB, but in the septum the normal value of Wext was not reached, which explains deterioration of function and reduction in LV stroke work (Fig. 8). During acute CRT, Wext increased for the septum and decreased for the LVfw. A nearly normal distribution of Wext was obtained during sustained CRT.
MEC for epicardial, midwall, and endocardial cell types.
To investigate the sensitivity of MEC in our model for variation in specific conductance of potassium currents, we performed an additional series of simulations with MEC in which Gto and GKs were set to values for epicardial, midwall, and endocardial cell types as proposed by ten Tusscher et al. (55). The only difference between epicardial and endocardial cell types is the value for Gto. As can be seen in Fig. 2, the smaller amount of Ito in the endocardial cell type leads to a smaller peak ICaL, but there is no significant difference in Ca2+ transient and Fnorm. Also with MEC, no significant differences in Ca2+ transient, Fnorm, or cardiac function between endocardial and epicardial cell types were observed (not shown).
Figure 10 shows Vmem, [Ca2+]i, and Fnorm for epicardial and midwall cell types. The lower value for GKs in the midwall cell type results in an increased calcium transient and a larger difference in contractile force between the septum and the LVfw during acute LBBB. During sustained LBBB, contractile force generated by the LVfw is still larger than contractile force generated by the septum, indicating that MEC has a somewhat different effect compared with the epicardial cell type. Regarding cardiac function, for both cell types the same trend as described for our setting of GKs was observed: MEC lead to a decrease in LV stroke work and an increase in EDV compared with acute LBBB. During sustained CRT, LV stroke work increased, and EDV decreased compared with acute CRT (not shown). Thus, regardless of the value for GKs, LV function deteriorated during sustained LBBB, whereas it improved during sustained CRT.
Our simulations illustrate that long-term cardiac MEC may explain the gradual worsening of cardiac function during LBBB as well as the gradual improvement of cardiac pump function after onset of CRT, whereas our predicted molecular and electrophysiological changes are in line with experimental and clinical observations. Therefore, MEC may serve as a framework for a more integrative understanding of remodeling processes involved in the asynchronous and the resynchronized heart.
In dyssynchronous and resynchronized hearts, MEC can explain a number of changes, such as the regionally different changes in APD, the change in dispersion of repolarization, and the change in pump function (stroke work, EDV) during chronic LBBB and CRT. In canine models of sustained LBBB, worsening of cardiac pump function as well as APD changes have been observed (51, 59). Regionally different changes in ICaL were seen in dogs with tachypacing-induced heart failure in combination with LBBB (2). There are also experimental data on the effects of CRT, including acutely improving pump function followed by ongoing improvement in the long run (58). Yu et al. (61) showed in patients with CRT a continuing increase in pump function following the acute hemodynamic effect of CRT. They also found a two-stage decrease after terminating CRT 3 mo after its start. Also, Delnoy et al. (14) showed gradual increase of ejection fraction during a 2-yr follow-up study of CRT. Aiba et al. (2) showed in their dyssynchronously failing hearts that CRT normalized the distribution of ICaL, calcium transient, and APD. Thus a large number of model predictions are supported by experimental and clinical evidence, supporting the validity of our model.
In a previous study, we used an earlier version of our model representing only the left heart and systemic circulation but with the same form of MEC. With that model, we were able to explain, not only the relatively uniform distribution of strains despite the moderate extent of asynchrony in electrical activation, but also the concordant T-wave during normal sinus rhythm as well as the T-wave changes during T-wave memory (20). In the present study, we found that MEC also leads to less dispersion in repolarization and better cardiac function in sustained CRT. Although an electrocardiogram cannot be deferred from this simplified model, increase in APD-60mV in early-activated regions and a decrease in APD-60mV in later-activated regions as shown in Fig. 6 would likely lead to a less concordant (or possibly discordant) T-wave. The predictions on electrical remodeling in CRT are supported by recent clinical studies. In one study, it was found that the amplitude of the T-wave decreases with longer duration of LBBB in patients (48), whereas our group confirmed electrical remodeling in patients with CRT by demonstrating the occurrence of T-wave memory in these patients (60).
Possibly the most interesting and intriguing result of the present simulation study is that a local electromechanical adaptation process, designed to create uniform distribution of mechanical myofiber workload, can lead to worsening of global pump function in the asynchronous heart. The explanation for this paradoxical finding is that, during the severe dyssynchrony in LBBB, the MEC-induced reduction in external work in late-activated LV free wall regions is not fully compensated by the increase in work in the opposing early-activated septum. This property is related to the assumption that GCaL can only change to a moderate degree in response to altered mechanical work. This assumption is based on experimental observations (2) that, compared with control dogs, ICaL in canine hearts with DHF was ∼20% higher in the early-activated anterior wall and ∼20% lower in the late-activated lateral wall. Also, our previous modeling study (20) revealed that allowing a maximum deviation of 25–50% around the normal value leads to estimations of electrophysiological and functional parameters that are close to clinical observations (20). Although in our previous study it was shown that a full compensation of pump function can only be achieved by a multifold increase in GCaL in early-activated regions, limiting the maximal adaptation to ±20% as used in the present study, will alter the outcome in pump function quantitatively but not qualitatively.
Closer observation of the model predictions shows why an apparently functional system of adaptations, like creating uniform distribution of external work, has adverse effects in LBBB hearts. In these hearts, MEC causes a reduction in ICaL in the late-activated, hyperdynamic LVfw. The septum, being approximately half the size of the LVfw, cannot compensate for the loss of work in the LVfw (Fig. 6) because this would require a much larger increase in Wext than can be achieved by an increase in ICaL within physiological range (2, 20). Experimental studies demonstrated that sustained LBBB leads to hypertrophy of the late-activated lateral wall (58, 59). Although we did not take this kind of structural adaptation into account in our simulations, it should be noted that a change in muscle mass alone cannot explain the electrophysiological changes predicted by the model and observed in the animal and clinical studies as described above.
The model simulations also demonstrate that there may be two mechanisms of action for CRT, acute recoordination of contraction and long-term electromechanical remodeling. After all, acute CRT improved cardiac pump function, whereas contraction became more disperse, and further improvement occurred without further resynchronization. Given the observation that the chronic effect of CRT may be as large as the acute effect, at least regarding ejection fraction (14), this study indicates the significant importance to further investigate electromechanical remodeling in CRT hearts. Summarizing, it appears that MEC can explain various data on electrophysiology, contraction, and pump function in dyssynchronous and resynchronized hearts.
Considerations regarding the modeling approach.
To the best of our knowledge, our model is the first to integrate cardiac electrophysiology, mechanics, and hemodynamics as well as MEC in a biventricular cardiac model. The model of the human action potential used in the present simulation study is the 2006 model of ten Tusscher and Panfilov (TP2006) (55), which is often applied in large-scale models of human cardiac electrophysiology. Other models of the human ventricular action potential are Priebe-Beuckelmann (PB1998) (41), ten Tusscher-Noble-Noble-Panfilov (TNNP2004) (54), Iyer-Mazhari-Winslow (IMW2004) (23), Bueno-Orovio-Cherry-Fenton (BCF2008) (11), Grandi-Pasqualini-Bers (GPB2010) (18), and O'Hara-Virág-Varró-Rudy (OVVR2011) (38). Elshriff and Cherry (16) compared action potentials of epicardial formulations of all these models. They concluded that all models differ in formulation and are usually based on different experimental data sets from humans. For the present simulation study, a good representation of Ca2+ dynamics is important, as it is influenced by changes in the L-type Ca2+ current. TP2006, GPB2010, and OVVR2011 describe ion concentrations and state variables for different subspaces. TP2006 contains a description of calcium dynamics in the cytoplasmic subspace and sarcoplasmic reticulum. Calcium in the subspace is buffered, and a diffusion flux is included to allow calcium released in the subspace to flow to the bulk cytoplasm. L-type Ca2+ current and ryanodine calcium release current inject calcium into the subsarcolemmal subspace, and, in turn, their dynamics are influenced by the subspace calcium concentration (55). GPB2010 and OVVR2011 have similar descriptions for calcium handling. Action potential shape and duration of the TP2006, GPB2010, and OVVR2011 models are comparable, but Ca2+ transients differ significantly among the three models. Peak Ca2+ concentration of the GPB2010 model is 0.4 μM, whereas it is about 0.8 μM for OVVR2011 and TP2006 (epicardial cell type).
Because of the close coupling between Ca2+ concentration and generated force, the mechanical behavior of our model may be significantly influenced by choosing a different model. Of the three models, TP2006 and OVVR2011 produce peak Ca2+ concentrations that lead to physiological contractile forces. Both in TP2006 and in OVVR2011, inactivation of L-type Ca2+ current is related to the calcium concentration in the cytoplasmic subspace. Because Ca2+ transients of the TP2006 and OVVR2011 models are similar, it is expected that mechanical behavior of our model would not significantly change if the TP2006 model would be replaced with the OVVR2011 model. Thus, for the present simulation study, similar effects of MEC with respect to cardiac function would have been found if OVVR2011 would have been used instead of TP2006.
Contraction of the sarcomere was modeled by model 5 of Rice et al. (46). Like other sarcomere models, this model is based on experimental data of rat trabeculae recorded at room temperature. Time-dependent contractile behavior of sarcomere models depends on, not only the amplitude and shape of the Ca2+ transient, but also on passive stiffness (collagen and titin), sarcomere length, and sarcomere shortening velocity. These influences are taken into account in our model. The calcium transient is obtained from the 2006 model of ten Tusscher et al. (55). A similar approach was taken by Campbell et al. (12), who integrated the model of Flaim et al. (17) of canine ionic currents and calcium handling with the 2008 model of Rice et al. (45). Unfortunately, models describing contractile behavior of myocytes from larger species are at present not available. Because in our model the L-type current is affected as a result of MEC, which in turn influences the calcium transient, changes in mechanical and hemodynamical behavior in our model are to a large extent related to changes in the calcium handling. Because of a lack of experimental data, it is at present not possible to validate the integrated cell model against human data. However, as discussed above, cardiac mechanics and hemodynamics correspond well with clinical measurements.
As one of few models, our model describes the entire systemic and pulmonary circulation, providing realistic boundary conditions for cardiac mechanics that are of particular importance for the simulation of long-term MEC. Compared with other whole-heart models, our model lacks detail in ventricular geometry. However, anatomical details are assumed to be of minor importance for a proof-of-principle study, as is the present study. The advantage of the relatively small number of segments (570 in total) was that the simulation of 1,000 heart beats to simulate both (sustained) LBBB and (sustained) CRT took less than 10 h on a high-end PC. A side effect of our modeling approach is the nonrealistically sharp changes in APD-60mV for the same values of GCaL observed in Fig. 6 at the junctions of RV free wall, septum, and LV free wall.
The idea that mechanically induced electrical adaptations may be involved in the regulation of cardiac electrophysiology is not novel. However, whereas many investigators use the term MEC for processes that occur within a single heart beat or several beats at most, the kind of MEC proposed in our study is a much slower process, presumably acting in the order of hours to days. This kind of MEC may, for instance, underlie the differences in vector ECGs between short-term and long-term LBBB observed by Shvilkin et al. (48).
In the present study, it is assumed that adaptation of GCaL is triggered by changes in (external) workload. Indeed, experimental observations suggest that changes in workload induce electrical remodeling (24, 50), but the exact trigger is not yet clear. Importantly, the present study does not claim any specific sensing mechanism. Local stretch, rather than workload, might be the trigger for MEC. After all, the amount of early-systolic prestretch is related to workload because it induces the Frank-Starling effect locally. However, prestretch by itself is not affected much by electrical remodeling (Fig. 5) and therefore is not likely to function as a regulating mechanism. Also, simulations using local stretch as trigger for MEC did not lead to stable solutions (E. Hermeling, T. Delhaas, F. Prinzen, and N. Kuijpers, unpublished data). On the other hand, Wext is a result of both stress and strain. It has been suggested that changes in stress and strain are sensed by Z-discs and influence cellular signaling molecules, including protein kinase C and A that regulate L-type calcium channels (43). An alternative feedback signal that can lead to similar results is local energy consumption. Recently, it has been shown that cardiac mitochondria may play a role in electrophysiology and calcium handling in response to oxidative stress (3).
Our choice for GCaL as a regulating factor for MEC does not exclude a possible role for other ionic membrane currents. In animal experiments using long-term asynchronous activation by ventricular pacing or LBBB, several changes in ionic membrane currents have been observed, including changes in ICaL (40), Ito (62), IKr, and IKs (37). In addition to changes in ICaL and calcium handling, Aiba et al. (2) observed differences in several potassium currents, including IK, IKl, and Ito in their DHF model. Although the differences in ICaL and calcium transient completely disappeared with 3 wk of CRT, changes in potassium currents were only partly restored. These results suggest that the effects of CRT may be more prominent in calcium handling than in potassium currents. In addition, Aiba et al. (2) reported that changes in ICaL and calcium handling were regionally different in DHF, whereas changes in potassium currents were not significantly different between regions. Our choice for modification of the specific conductance of ICaL in the present simulation study was based on these experimental findings. The relatively large effect of small changes in ICaL on the calcium transient and contractile force in our model is explained by the increase and decrease in storage of calcium in the sarcoplasmic reticulum during diastole, which, in turn, affects calcium release during systole. This secondary effect of changes in ICaL accumulates over several heart beats.
Many changes occur in dyssynchronous hearts, among others in sarco(endo)plasmic Ca2+-ATPase (SERCA), RyR, and the NCX (9). These changes influence calcium handling and thus contraction. However, Aiba et al. (2) did not find significant regional differences in SERCA, RyR, and NCX. Because our concept of MEC requires regionally different responses, we did not include changes to SERCA, RyR, and NCX in the present simulation study. In addition, metabolic dysfunction has been reported by Agnetti et al. (1) in the late-activated LV lateral wall, but, to our knowledge, no data are available on regional differences. For this reason, metabolic dysfunction was not considered in the present simulation study.
In previous simulation studies, we have shown that inclusion of transmural heterogeneity of Ito and IKs (25) or remodeling of IKs instead of ICaL (20) also leads to a uniform distribution of external work and T-wave concordance. Also in the present study, we performed a series of simulations in which GKs was adapted instead of GCaL. In these simulations, GKs could vary in between 0.098 nS/pF and 0.392 nS/pF. When GKs was adapted, changes in action potential duration were much more prominent than when GCaL was adapted (Fig. 11), but the effect on cardiac function was similar (not shown). Plotnikov et al. (40) observed a shift in ICaL activation and slower inactivation in canine epicardial myocytes after 3 wk of pacing from the LV epicardium. The assumed long-term MEC in this study does not specify the exact mechanism of the change in ICaL, i.e., number of channels, phosphorylation, or channel kinetics.
In our model, it is assumed that all segments within one ventricular wall bear the same amount of stress. This approach is valid when the ventricular walls do not contain irregularities in mechanical properties such as infarcted regions or scar tissue (4, 5). With this relatively simple geometry, we have recently shown that strain patterns can be obtained that agree with clinical measurements of strain patterns in patients with LBBB (32).
Simulations with our model of the heart and circulation predict that electromechanical remodeling reduces the amount of work generated by the late-activated LVfw, but this decrease is not fully compensated by an increase in Wext generated by the early-activated septum. This results in the paradoxical situation that, during longer-lasting dyssynchrony, electromechanical remodeling causes a decrease in pump function despite the reduction in dispersion of mechanical workload and repolarization. Conversely, electromechanical remodeling optimizes cardiac pump function in less dyssynchronous ventricles. Consequently, the model pinpoints two complementary mechanisms of improved pump function by CRT, acute recoordination of contraction and long-term electromechanical remodeling.
This work was partially funded by a personal research grant within the framework of the Dr. E. Dekker program of the Dutch Heart Foundation (NHS-2012T010 to J. Lumens).
No conflicts of interest, financial or otherwise, are declared by the authors.
Author contributions: N.H.K., E.H., and F.W.P. conception and design of research; N.H.K. and E.H. performed experiments; N.H.K., E.H., J.L., H.t.E., T.D., and F.W.P. analyzed data; N.H.K., E.H., J.L., H.t.E., T.D., and F.W.P. interpreted results of experiments; N.H.K., J.L., and T.D. prepared figures; N.H.K. and F.W.P. drafted manuscript; N.H.K., E.H., J.L., H.t.E., T.D., and F.W.P. edited and revised manuscript; N.H.K., E.H., J.L., H.t.E., T.D., and F.W.P. approved final version of manuscript.
Appendix A: COMPUTATIONAL ASPECTS OF ELECTROPHYSIOLOGY
The bidomain model can be written as a coupled system of ordinary differential equations and partial differential equations (27). To compute Vmem, the differential equations were solved using a forward Euler scheme. To compute Vext, a system of linear equations was solved each time step by an iterative method as described in Kuijpers et al. (27). Criteria for segment size were obtained by cable theory and considering subthreshold behavior along the cable as previously described (27). For the bidomain parameters in Table 1, a length constant of 0.09 cm was found. For proper simulation of action potential propagation, segment size was set to 0.01 cm. Note that this segment size was chosen for accurate simulation of impulse propagation and is not related to the physical dimensions of the ventricular walls.
Because we used an explicit forward Euler scheme, simulation time step Δt had to satisfy the stability condition as formulated by Puwal and Roth (42). For a serial arrangement of electrically coupled segments, this leads to (11) where Δx is segment size (0.01 cm). For the bidomain parameters defined in Table 1, this condition is satisfied when Δt < 0.13 ms. In the present study, we used Δt = 0.01 ms.
Ionic currents and concentrations were computed using a forward Euler scheme with a time step of 0.01 ms. Gating variables for the ionic currents were computed using the Rush-Larsen method (47). Computation time was reduced by increasing the simulation time step for the ionic currents and concentrations to 0.1 ms during repolarization and rest. This led to a reduction of 70% in computation time, without significant loss of accuracy (30).
A well-known deficit of membrane models is drift in intracellular Na+ concentration ([Na+]i) and intracellular K+ concentration ([K+]i) during longer-lasting simulation runs. In the present study, simulation runs of up to 1,000 cardiac cycles were performed. To prevent [Na+]i and [K+]i from drifting, a regulating feedback step was performed each simulation time step Δt as proposed by Van Oosterom and Jacquemet (56): (12) (13) with [Na+]i,rest = 8.3 mM and [K+]i,rest = 137.25 mM (55).
Appendix B: COMPUTATIONAL SCHEME FOR MYOFIBER MECHANICS
Ca2+-force relations as well as cardiac mechanics and hemodynamics were computed using a forward Euler scheme with a time step of 0.01 ms. Ventricular mechanics and hemodynamics in our model were solved as described by Lumens et al. (33). To solve cavity mechanics, a computational scheme for a serial arrangement of mechanically coupled segments was required that allowed applying changes in overall strain with and without the advancing of simulation time. We recently published a new computational scheme for fiber mechanics (20) that was based on our previously published scheme (29). This scheme was further refined to compute myofiber mechanics for the biventricular model of the present study.
The numerical scheme to solve the equations for the three-element mechanical model was based on the scheme introduced by Solovyova et al. (49). The mechanical state of each segment was defined by lsi, ls, and Fnorm. Each simulation time step, first lsi was updated using velocity of internal sarcomere shortening v = −dlsi/dt. Next, Fnorm was computed using lsi as well as [Ca2+]i obtained from the model of ten Tusscher et al. (55). Finally, ls was computed by solving Δls as described below.
In this derivation, a fiber is defined as a serial arrangement of scaled segments. All segments have equal reference lengths l0 (expressed in cm) corresponding to reference sarcomere length ls0 (expressed in μm), i.e., l0 = ξ·ls0, where ξ is a scaling factor. First Piola-Kirchhoff stress Tsegment generated by a single segment is equal to first Piola-Kirchhoff stress Tfiber generated by the fiber. Thus also ΔTsegment is equal to ΔTfiber. Δls for each segment is found by solving a system of equations as follows. For each individual segment, it holds Tsegment = TSE + TPE, which implies (14) where lSE denotes the length of SE defined by lSE = ls − lsi and (15) (16) (17) (18)
Rewriting Eq. 14 and using ΔlSE = Δls − Δlsi = Δls + v·Δt gives an expression for Δls in ΔTsegment, Δt, and ΔFnorm: (19) where α, β, and γ are defined by (20) (21) (22)
Actual length of scaled segment n is denoted by ln and is defined by ln = ξ·ls,n. Assuming that a fiber is composed of N segments coupled in series, fiber stretch λfiber is defined by (23) where L and L0 denote the actual fiber length and reference fiber length (unit cm), respectively. From equation (23), it follows that (24)
From Eq. 25, ΔTsegment and ΔTfiber are computed by (28)
Internal sarcomere shortening velocity v needed to compute Δls from Eq. 19 is obtained as follows. Rewriting Eq. 3 and inserting definition 4 gives for each individual segment (29) from which v can be obtained by (30)
Here, TCE is computed from Eq. 5 and 8 by (31)
Initially, it is assumed that the electrophysiological state of all segments is resting and that Fnorm = 0. Thus TSE = TCE = 0 and lsi = ls for all segments. Stress generated by the segment must come from the parallel elastic element and is related to λfiber = exp(ϵfiber), with ϵfiber natural strain corresponding to the initial cavity volume.
Each simulation time step, cavity volumes are updated as follows: 1) new cavity volumes for RV and LV are obtained by integrating blood flow through the valves [see Lumens et al. (33)]. 2) Natural fiber strain (ϵfiber) for RVfw, septum, and LVfw is obtained from the new RV and LV cavity volumes by applying an iterative process with Δt = 0 to reach equilibrium of forces at the junction of the three wall segments [see Lumens et al. (33)]. In each iteration, the computational scheme described below was applied with Δt = 0 to compute the sarcomere length ls for all segments. 3) Finally, the computational scheme described below was applied once more with Δt, the actual simulation time step.
To update the mechanical state of all segments for the newly found natural fiber strain ϵfiber, the following scheme was used, either with Δt = 0 (step 2 above) or Δt > 0 (step 3 above): 1) Δλfiber = exp(ϵfiber) − λfiber is computed from the new fiber strain ϵfiber and λfiber from the previous iteration. 2) For each segment, internal sarcomere length lsi is computed using Δlsi = −v·Δt, with internal sarcomere shortening velocity v computed during the previous iteration. Note that Δlsi = 0 when Δt = 0. 3) For each segment, Fnorm is updated using lsi and [Ca2+]i obtained from the model of ten Tusscher et al. (55). This is only done in case Δt > 0. 4) αfiber, βfiber, and γfiber are computed from Eq. 26–28 using αn, βn, and γn from the previous iteration, and Fnorm. 5) First Piola-Kirchhoff stress ΔTfiber and ΔTsegment are computed from Eq. 29. 6) For each segment, sarcomere length ls is computed from Eq. 19. 7) For each segment, first Piola Kirchhoff stress for SE and PE (TSE and TPE) are computed from Eq. 5–7. 8) For each segment, v, α, β, and γ are computed from Eq. 31 and Eq. 20–22. 9) Finally, λfiber is computed from Eq. 23.
- Copyright © 2014 the American Physiological Society