Myocardial infarction (MI) significantly alters the structure and function of the heart. As abnormal strain may drive heart failure and the generation of arrhythmias, we used computational methods to simulate a left ventricle with an MI over the course of a heartbeat to investigate strains and their potential implications to electrophysiology. We created a fully coupled finite element model of myocardial electromechanics consisting of a cellular physiological model, a bidomain electrical diffusion solver, and a nonlinear mechanics solver. A geometric mesh built from magnetic resonance imaging (MRI) measurements of an ovine left ventricle suffering from a surgically induced anteroapical infarct was used in the model, cycled through the cardiac loop of inflation, isovolumic contraction, ejection, and isovolumic relaxation. Stretch-activated currents were added as a mechanism of mechanoelectric feedback. Elevated fiber and cross fiber strains were observed in the area immediately adjacent to the aneurysm throughout the cardiac cycle, with a more dramatic increase in cross fiber strain than fiber strain. Stretch-activated channels decreased action potential (AP) dispersion in the remote myocardium while increasing it in the border zone. Decreases in electrical connectivity dramatically increased the changes in AP dispersion. The role of cross fiber strain in MI-injured hearts should be investigated more closely, since results indicate that these are more highly elevated than fiber strain in the border of the infarct. Decreases in connectivity may play an important role in the development of altered electrophysiology in the high-stretch regions of the heart.
- left ventricle
- mechanoelectric feedback
damage resulting from myocardial infarctions (MIs) can cause changes in the overall structure of the heart and significantly alter its function. In many cases, such damage leads to a cycle of progressive change that continues to reduce the efficiency of the heart, ultimately ending in heart failure (HF) (32). In addition, the structural changes to the myocardium after an MI can also alter the electrical activity of the heart, potentially leading to dangerous arrhythmias that, at their worst, end in sudden cardiac death (SCD) (29). Because MIs are a leading cause of both HF and SCD, and these two conditions create an immense burden both socially and economically (1), understanding how an MI creates the changes in ventricular structure and function that lead to these clinical endpoints is an important target for continued research and development.
Mechanical function and electrical activity are closely linked in the heart, and these two factors play an important role in the postinfarct heart. Specifically, changes in regional structure and material properties after MI (24) alter regional strain and stress (40) and can decrease pump function. Altered patterns of stretch and load may be causal in the continued remodeling of the ventricle and ultimately lead to HF (10, 18). Meanwhile, as strain dynamically alters cellular currents that drive electrophysiology through mechanoelectric feedback (MEF) (8, 15), changes in mechanics that result from MIs may be a driver in creating electrical instabilities that can lead to SCD.
Although experimental and clinical studies of electromechanics in the heart have advanced greatly, allowing detailed measurements of mechanical and electrical activity of the myocardial tissue, notable gaps continue to exist in what can be measured. For example, while magnetic resonance imaging (MRI) (31) and echocardiography (35) can provide detailed strain analysis of the myocardial wall, due to the complex anisotropy and architecture of the heart and the nonlinear stress-strain relationship of myocardial tissue, directly measuring the mechanics of the oriented sheets of muscle cells is not possible. In addition, optical mapping provides relatively easily obtainable direct experimental measurements of electrical activity on the surfaces of the heart (9), yet measurements of electrophysiology inside the heart wall are very difficult to obtain experimentally.
One technique that can help fill these gaps in experimental research is the use of computational methods. Computer simulations have proven themselves to be a useful tool in understanding both mechanical and electrical activity of the heart, and, quite recently, models have advanced to sufficient complexity to tackle coupled problems in electromechanics (2, 14, 22). The use of computer simulations provides theoretical insight into heart function that can be validated later through experimental work, or can provide a means for testing novel hypotheses in function and/or clinical treatment (41, 46). The purpose of this study was to theoretically examine the myofiber strains in a left ventricle (LV) with an MI using computer simulations and detailed in vivo measurements. Using a coupled model of electromechanics, we were able to evaluate the entire cycle of the heart, highlighting regions of time that are of mechanical interest. The myocardial strains were then examined for the role in changing the electrical dynamics of the heart through MEF. In addition, the connectivity between cells in the border zone was altered as well to examine how known processes such as fibrosis could exacerbate changes to electrophysiology.
Mathematical model of the infarct injured LV.
A mathematical model describing the myocardium was used to simulate the function of the LV with an MI through multiple cardiac cycles. This model consists of a set of tightly coupled subsystems expressing continuum mechanics, electrical diffusion, and cellular dynamics, which describe how changes at the cellular level give rise to active force generation and how the motion of the heart changes as a result of hemodynamic forces and active contraction (38).
Briefly, tissue mechanics is modeled using the equations of equilibrium formulated with respect to the undeformed geometry in terms of the deformation gradient (F) and the second Piola-Kirchhoff (2nd PK) stress tensor (S). This stress is divided into a passive term (Sp), and an active component (Sa). The passive component is calculated through a strain energy function (Ψ) and the Green Lagrange strain (E) describing a nearly incompressible, Fung type transversely isotropic hyperelastic material (11, 12). The active 2nd PK stress is converted from a Cauchy stress (σa), which is calculated from a reference tension (To), the primary fiber stretch ratio (λ), the fiber stretch velocity (), the current cellular state (s), and the Jacobian (J). 1 2 3 4 The cellular state, s, is a vector holding values related to the ion cycling and force generation in myocytes. It is described using a series of 33 ordinary differential euqations (ODEs) (16, 43) describing transmembrane and intracellular ion currents, calcium binding to troponin C, and subsequent exposure of cross-bridge binding sites: 5 where s holds state variables at time t, υ is the transmembrane voltage, and λ is the primary fiber stretch ratio.
The spatiotemporal dynamics of the voltage (υ) through the myocardium is determined by a bidomain model that describes ion flow through the intracellular and extracellular myocardial spaces (27), 6 7 Here Iion is the total ionic current scaled to membrane capacitance, Mi and Me are the intracellular and extracellular conductivity tensors, respectively, and ue is the extracellular potential.
This model is solved using the finite element method in Cartesian coordinates using standard methods (25). Operator splitting is used to solve the coupled problem, where at each time step 1) the ODE system is solved updating the state vector s, 2) the bidomain partial differential equation is solved updating the resulting voltage field, and 3) the mechanics is solved using the most recent values of the state vector s.
Geometric model of the ovine LV with MI.
A geometric mesh representing the ovine LV was constructed from previously obtained MRI data (40) from a sheep that had undergone surgical ligation of the left anterior descending coronary artery to establish a fully dyskinetic anteroapical infarct. This procedure had created a transmural aneurysm of ∼30% of the ventricle wall at 22 wk postprocedure, and the resulting heart suffered from a severely depressed ejection fraction. Long- and short-axis tagged MRI images were used to contour the endocardium and epicardium at minimum volume during the cycle (Fig. 1A), and these contours were used to construct a geometric mesh of 18 × 12 × 1 (longitudinal × circumferential × transmural) triquadratic 27-node elements with a small 9-element apical cap. This geometry was broken down into three distinct sections as previously described (13, 40): a normal remote myocardium, a hypocontractile border zone, and a noncontractile aneurysm (Fig. 1B). Fiber angles in the mesh were assigned from previous diffusion tensor imaging of the explanted heart (39), assigned at nodes and interpolated through the elements. Fiber angles were set to the circumferential direction in the infarct (Fig. 1C).
Ventricular mechanics were analyzed on this 2,775-node mesh. Passive material properties were chosen iteratively, starting from initial values published using a similar mesh constructed in prolate spherical coordinates (40). Subsequently, the stiffness scaling factor was increased to match the in vivo end-diastolic volume of the LV. In addition, active stress in the border zone was reduced 50% by changing the To parameter, in line with previous estimates made (13). The infarct region was set much stiffer than the functional myocardium and had active contraction set to zero.
The mechanics mesh was refined two times into an ∼130,000-node mesh on which the cell model and the diffusion of electrical signal were simulated. Cell models were not run in the infarct, and electrical diffusion was reduced by a factor of 1,000 in this acellular region, confining the cellular activity and electrical signal propagation to the border zone and healthy regions.
Cardiac cycle and model boundary conditions.
The cardiac cycle was modeled using a dynamic endocardial pressure and an activating current against an assumed constant epicardial pressure. Starting from an unloaded ventricle with cell nodes at the polarized state of −95 mV, the pressure was increased on the endocardial surface until it reached the in vivo end-diastolic pressure that had been measured in the animal by an inserted catheter. Twenty milliseconds before reaching end-diastolic pressure, activation of the ventricle was simulated with a current added sequentially over 10 ms from the apex up along the lower two-thirds of the endocardial surface as a rough approximation of the Purkinje network. The effect of the infarct was ignored for activation, and, since the cellular model was turned off in the infarct, this region did not contribute to the activation pattern. The chosen activation sequence allowed propagation of the signal from apex to base and from endocardium to epicardium. After stimulation, as active force increased, the endocardial pressure in the ventricle was increased to maintain constant volume until the prescribed systolic pressure of 11.76 kPa was reached. The ventricle then contracted against a two-element Windkessel model until the aortic pressure exceeded the calculated ventricular pressure. Two-element Windkessel model parameters were taken from published ovine data (33) and modified within the given ranges until the model gave the correct end-systolic volume. Decreases in pressure with constant volume then progressed in isovolumic relaxation until t = 800 ms. This cycle was then repeated until a steady-state pressure-volume (PV) loop was reached. All simulations were run either locally on a laptop using a 2.4-GHz Intel Core 2 Duo central processing unit or on single nodes of a computational cluster using Intel Xeon 2.5-GHz processors. Each complete 800-ms loop took ∼20 h to solve on both systems.
We introduced MEF in the model with nonspecific cation stretch-activated currents (SAC). A linear SAC was added to the cell model (Eq. 7), having the form: 8 9 10 Here msac is the magnitude of the slope with respect to the current, gsac is the maximum conductance, and ER is the reversal voltage. In our model, although there is a range of reversal potentials reported between 0 and −11 mV (4, 21, 45), we chose a constant reversal voltage of −6 mV as reported in rat ventricular myocytes (45). Because the role of stretch-activated channels and mechanoselective gated channels in the myocardium is still very unclear, we ran a number of simulations to evaluate the trend of increasing current with increasing stretch. In our model, we did this by varying the linear response (gsac) of this current in a series of simulations from 0 to 250 (S/Fs), which corresponds up to two times the published values for these types of currents (15, 45), to account for other possible similar current sources and establish potential trends.
In addition, simulations were run where the electrical conduction in the border zone was scaled to examine the role of stretch-regulated currents in relation to dysfunctional electrical connectivity. Full-cycle simulations were run with bidomain intracellular and extracellular connectivity reduced by a factor of 10×, 100×, and 250× to simulate potential fibrosis or other breaks in cellular connectivity due to remodeling after the infarct.
Base model of the injured ventricle and comparison with in vivo measures.
The coupled electromechanical simulation captured all four phases of the cardiac cycle and reached steady action potentials, calcium transients, pressure histories, and PV loops after three cycles (Fig. 2, A–D). End-diastolic volume and end-systolic volume of the final steady state were matched to in vivo measurements at the prescribed pressures. In addition, after the simulation reached a steady state, the resulting minimum ventricular pressure of ∼2 kPa was close to the published minimum of 3.58 kPa that had been directly measured in this animal (40). Simulated ejection was somewhat faster than experimental measures, with end ejection occurring at ∼460 ms into the simulated cycle compared with ∼520 ms as observed with MRI.
Global cartesian strain in the myocardium was calculated through systole on the mechanics mesh, and circumferential strain in the remote region and border zone during systole was consistent, but the contraction was slightly faster than MRI measurements (Fig. 3, A and B). At end systole, the global circumferential strain was consistent with MRI measurements (Fig. 4, A and B) with very heterogeneous strain around the infarct. Individual strain components were also compared, and, in the region of the border zone, especially good agreement of the normal circumferential and radial strain components was observed in cross section (Fig. 4, C and D). Discrepancy was seen in the longitudinal direction (Fig. 4E) and shear components that used this orientation (Fig. 4, F–H). However, this strain component has known limitations when measured using MRI (5, 6), and it is unknown whether this discrepancy is due to the model or the measurement technique.
Mechanical results of the LV with an MI.
Calculated fiber and cross fiber strain cycles for the steady-state PV loop are depicted for the remote, border zone, and infarct regions (Fig. 5). At each element, the midwall strain was calculated, and the elements grouped into the three regions. Average of strain in each of the groups and the SD is shown for these regions over the time course of the beat. As expected, most elements underwent extension in the fiber and cross fiber direction during diastole, with peak average end-diastolic strain of ∼0.05 for both strain components. After the initiation of activation, elements in the remote region begin to contract while the infarct continues to extend as pressure is equalized over the contracting ventricle and noncontractile infarct. The hypocontractile border zone between those regions behaved in a much more dispersed fashion, since areas closest to the infarct in the thinned region remain in extension while those in other areas begin to contract rapidly. During ejection, this trend continues, as remote regions continue to contract while infarct and near-infarct border zone regions stay in extension until the pressure ramps down in isovolumic relaxation. The cross fiber strain (Fig. 5, D–F) is significantly more dispersed in behavior than the fiber strain (Fig. 5, A–C). While a fraction of the near-infarct border zone regions remains in fiber extension for the entire cardiac cycle, in cross fiber strains, there are elements that rapidly expand during systole, and the average strain remains in extension for the entire border zone.
Electrical dysfunction of the LV with an MI.
Figure 6 depicts the effect on the addition of stretch-activated channels in the absence and presence of border zone connectivity reduction. Two different areas in the mesh are tracked, one in the remote region and one in the border zone, with differing fiber strain profiles over the cardiac cycle (Fig. 6A). With the addition of SAC currents and their increasing magnitude, there was a very slight decrease in observed action potential duration in the remote region, on average decreasing ∼7 ms from the lowest to highest SAC magnitude (Figs. 6B and 7B). Meanwhile, there was a corresponding decrease in the spread of the action potentials, with the SD decreasing from ∼13 ms to ∼5 ms with increasing SAC magnitude (Fig. 7D).
In the border zone region, however, action potentials in regions of continued strain increased slightly with increasing magnitude of SACs. When decreased electrical connectivity was also included in the model, the observed trends in the border zone region increased. With electrical decoupling of nearby nodes, the action potential lengthened substantially in areas of continuous stretch (Fig. 7A), increasing in average from ∼322 to 330 at the highest SAC strength. In addition, the dispersion of the action potential duration increased dramatically, going from below 5 ms in the control simulations to over 23 ms with the most dysfunction of the electrical diffusion and the highest-magnitude SACs (Fig. 7C). The response was highly nonlinear, with only small changes in action potential duration with increasing gsac until a large step observed at the highest magnitude (250 S/Fs), with dispersion of duration also stepping up dramatically at between 100 and 250 S/Fs.
These results provide, to our knowledge, the first strongly coupled electromechanical model of the LV with an MI. Our simulation has attempted to match measurements of volume, pressure, and strain from the beating ventricle and to provide useful insight into the function and dysfunction of the injured heart. The simulations allowed us to compute fiber strains through the entire cardiac cycle and to theoretically test their potential effects on electrophysiology. Our main results with this model are the calculation of an increase in fiber and cross fiber strains in the sensitive border zone during systole, with some tissue near the infarct remaining in fiber and cross fiber extension through the entire cardiac cycle. We also demonstrate that the heterogeneity of mechanical dysfunction in the border zone region can be especially detrimental to electrical activity in the presence of altered connectivity of the local myocardium. We found that, even with large increases in stretch-induced currents, changes to action potential are largely smoothed out in the myocardium under normal conduction velocities and that changes to the underlying cellular and extracellular connectivities were required for the development of abnormal action potentials.
We were able to use our model to broadly capture global mechanics, with good correspondence of simulated strains to those measured by MRI. The rate of contraction was slightly faster in the model than in the animal, but strains were largely consistent at end-systole. Individual strain components were also in good agreement, with the exception of the normal longitudinal strain component (ELL) and the shear circumferential longitudinal strain (ECL) and shear radial longitudinal strain components. Both the model and MRI show circumferential strain that varies from stretch to compression at end-systole based on the location of the infarct. Radial strain follows the same pattern. For ELL, however, the MRI data show the entire cross section in longitudinal compression at end-systole regardless of the circumferential location. This is unexpected, since a slightly incompressible material like myocardium should have extension in one normal component if the other two are in compression and points to possible error. As error in the MRI strain mapping in the longitudinal orientation is known to be the greatest, this measurement is potentially why the discrepancy between simulation and MRI is observed in the shear ECL component, where large phasic differences are observed between the simulation and the calculated MRI strain.
Despite these uncertainties, model agreement is sufficient to broadly examine fiber dynamics in our simulation. As has been reported previously (13, 40), we found that, in the ventricle with an MI, fiber strains are highly heterogeneous, which is expected in light of the experimental measures of mechanical dyssynchonry in injured hearts (34). We calculated a highly similar pattern of extension and contraction in fiber and cross fiber directions in the healthy region remote from the infarct in the ventricle, but strain profiles changed substantially as one moved though the border zone and into the infarct. The noncontractile aneurysm extends during systole and tethered to the border zone; our results show that the border zone will experience heightened strains as a result. Previous experimental studies of the hearts with infarcts have indicated stretch in the border zone during end systole (28), and these have been demonstrated in theoretical models as well (13, 40). However, by using a full model of the heart cycle, we now show that the greatest strain differences appear at the end of isovolumic contraction. In particular, we see very high peaks in cross fiber strain at the beginning of ejection.
It is postulated that abnormal mechanical forces in the border zone may be a driving force for continued remodeling of the organ and the progression of HF (10, 17, 18). Our results give new insight into the strain history of muscle fibers during the heart beat and may possibly be useful for targeting therapies and explaining their mode of action. For example, passive constraint devices, such as the Acorn Corcap device, have shown to be effective in mitigating the progression of HF, and support against increased strain and stress is hypothesized to be a mechanism of action (30). Finite element method studies have demonstrated that end-systolic stress should be minimized with the use of a passive restraint (19). However, our results indicate that strain is most dispersed at the beginning of ejection, and measures of how such devices and techniques minimize abnormal strains at the beginning of ejection may be more important than their effect at end systole.
One potential confounding factor to our simulations is the role of the right ventricle (RV) in the developed stress pattern, especially since during systole the RV may alter strain and stress compared with an LV alone. Activation patterns and RV dynamics certainly have an effect on LV dynamics, since pacing studies have shown that how the stimulating current develops through the two ventricles can cause measurable changes to left ventricular strain (44). However, because our model attempts to globally match strain dynamics, the predicted overall stress is likely similar when an RV is included. In addition, in our special case of an anteroapical infarct, the area of greatest interest is in the border zone that is somewhat distant from the septal region where any effect of the RV will be largest.
Effect on cellular dynamics.
With even a simple single current model of MEF, we were able to simulate many of the features observed in experimental testing. Our model simulated increased action potential duration with increased stretch in the border zone, in line with a wide range of animal and experimental models where trends of increasing action potential duration with stretch were observed (7, 23, 26, 42). In addition, we observe a slight reduction in action potential duration dispersion with increasing SAC activity, which is in line with experimental measures that show that stretch-activated channels in healthy tissue decrease action potential duration dispersion (36). Because all baseline dispersion in the absence of SACs occurs from diffusion differences in the model, while it is expected that duration should decrease slightly on average with SAC activity, the decrease in dispersion is somewhat unexpected. How a nonspecific cation channel can reduce action potential duration spread caused by diffusion differences in the ventricle geometry remains unclear but in line with previous studies (36).
We predict that the same mechanism can cause the opposite to happen in dysfunctional border zone myocardium, increasing action potential dispersion in the region and electrical heterogeneity, potentially making the area more susceptible to arrhythmias. In our model, there appears to be some resilience to ion flow changes caused by SAC activity, with only slight changes to electrophysiology as currents increase. However, at some magnitude, the cell model starts to produce a strong change in action potential, with a large step in average duration and dispersion.
Implications of the role of connectivity.
Although stretch has been shown to increase action potential duration both in single cell studies (42) as well as whole tissue preparations (26), in our simulations, we only see significant changes to action potential duration and dispersion in the presence of decreased connectivity. This could be because of the magnitude of SAC currents used in our model, since, even though we used currents that exceeded published values, there are indications that MEF is exacerbated in the heart when an MI is present (23, 26). Although there have been previous models that have generated arrhythmias in cardiac tissue with the application of SAC currents, more than fivefold increases in published SACs were required (20). Another potential explanation is that tissue preparations of cells or tissue are subjected to more homogeneous strain than the infarct model we are examining and that, under similar strains, cells will work in concert and any changes to electrophysiology will be readily detectable. In the heterogeneous strain fields of the heart with a developed aneurysm though, very few myofibers share the same strain profile, and the electrical diffusivity of the myocardium is sufficient to smooth out the aberrant electrophysiology of isolated pockets of dysfunction.
Nevertheless, we were able to create significant changes to action potential duration and dispersion through changes in connectivity of the myocardium. These results have particular application to the injured heart, since fibrosis is a common phenomena associated with MIs and remodeling of the ventricle (37) and can alter the electrical connectivity of the region. This may be especially important in the progression of HF, since fibrosis can extend into healthy myocardium during the process and could potentially create areas of decreased electrical connectivity that would then be more susceptible to strain effects.
While capturing many aspects of the dynamics of the LV cycle, limitations in accurately simulating the entire heart beat obviously remain. One important limitation is the choice of stimulating current used in the simulation. Without complete data on the activation of the infarct injured heart, we assumed an endo-to-epicaridium, apex-to-base activation pattern. There is the possibility that the infarct has disrupted this pattern substantially, and different activation patterns could alter ventricular mechanics results and the electromechanical feedback studied. In addition, the 225-element mechanics model is just a basic approximation of the ventricle, with numerous factors that have been ignored that may alter the results. One of these factors is the spatial resolution, since, while our LV model is divided into the three different materials, the individual regions are homogeneous, and we ignore possible heterogeneities in the myocardium, such as perfusion heterogeneities from the infarct that could alter local myocardial properties. In addition, no valve action is directly modeled, the effect of the pericardium is ignored, and only a basic circulation model is used in our simulations, meaning that the pressure cycle is only a simplification of the actual pressure the LV experiences. Further work is needed to develop models of higher spatial resolution and better boundary conditions and that include valve dynamics and better circulation coupling to create more realistic models that may be of greater use to clinicians and researchers.
An additional limitation in the analysis of these studies is the nature of the simulated MEF channels. We use a simple linear SAC with one chosen reversal potential, and, while nonspecific channels have been found in numerous species (21), the effect of strain and the effect of blocking of these channels depend greatly on the species used (15). Furthermore, changes to electrophysiology can occur from sources other than nonspecific SACs, such as the Na+/Ca2+ exchanger (3), indicating that our model of MEF feedback is still incomplete. Finally, even with well-understood SACs, the effect of the individual strain components on current is completely unknown. Because we see very large changes in cross fiber strain in addition to those seen in fiber strain after an infarct, how strain can accurately be coupled to electrophysiology remains unclear. Until there is a clearer experimental picture on the effect of stretch on cell dynamics, including a complete description of the many potential cellular pathways, the differential effects of strain components, and the potential effect of strain rate, substantial limitations in coupled electromechanical simulations employing MEF will remain.
Conclusions and future directions.
We have presented a strongly coupled electromechanical model of the LV with an MI and used this model to derive theoretical fiber strain profiles as well as study the potential implications of elevated strains in destabilizing the electrical substrate of the heart. Such mathematical models can provide an additional means of insight into the function and dysfunction of the heart and may be useful to derive hypotheses for experimental study or for use in optimizing therapies intended to reduce mechanical or electrical problems in the myocardium. Future work in this area may be done to compare the infarcted chamber with the noninfarcted chamber, as well as to move to human models of the myocardium, the use of biventricular meshes and better circulation models, and better numerical methods for solving the complex mathematical systems that develop from describing the multiscale physical description of the working myocardium.
The research was supported by a Center of Excellence grant from the Research Council of Norway (RCN) to the Center for Biomedical Computing, Simula Research Laboratory, and by RCN Grants YFF-162730 and FRIBIO-205381. We would also like to acknowledge National Heart, Lung, and Blood Institute Grants 5R01HL-063348, 5R01HL-084431, 5R01HL-086400, and 5R01HL-077921, which provided the in vivo data for simulation input and model validation.
No conflicts of interest are declared by the authors.
- Copyright © 2012 the American Physiological Society