| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1Department of Biomedical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands; 2Laboratory of Cardiac Energetics, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, Maryland; and 3Department of Physiology, Maastricht University, 4Department of Lead Modeling, Medtronic Bakken Research Center, and 5Department of Biophysics, Maastricht University, Maastricht, The Netherlands
Submitted 5 April 2005 ; accepted in final form 12 June 2005
 |
ABSTRACT |
|---|
 TOP ABSTRACT MATERIALS AND METHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
eikonal-diffusion equation; electromechanics; ventricular pacing; finite elements
These effects are most likely due to abnormal and asynchronous electrical activation of the ventricles. When pacing from a site in the ventricles, depolarization spreads more slowly and less uniformly (21, 35, 38) than normal. As a result, the distribution of local amplitude and time course of contraction is affected considerably (22, 38) and was found to vary systematically with increasing distance from the pacing site (6, 19, 22). The goal of this study was to investigate whether the regional time courses of strain, as measured, may be sufficiently explained by coordinated actions of known mechanical, electrical, and structural properties of normal myocardial tissue. For that purpose we compared a simulation of cardiac contraction, based on the latter properties, with direct measurement of cardiac contraction, using MRI tagging in experiments on dogs. Thus we may obtain better insight into mechanisms relevant to pacing, with the final goal of improving ventricular pacing in general.
Several computational models of cardiac electromechanics have been used to investigate the effect of depolarization sequence on cardiac contraction globally in a normal heart (28), regionally in a normal heart (13, 30), during ventricular pacing (15, 32), and with regional ischemia (18). Vetter and McCulloch (37) investigated the influence of stretch on depolarization, but active contraction was not included. Rice et al. (25) and Sachse et al. (27) designed a more comprehensive model of cardiac electromechanics. These investigators coupled a canine (25) and a human (27) cellular ionic model to a model of myofilament activation (26). However, this was not incorporated in the whole heart because of computational demands. In the study of Usyk and McCulloch (31) the influence of biventricular pacing was investigated in a failing heart.
The model used in the present study was originally used to simulate a normal heartbeat (13). Next, we estimated the model's parameter values for electrical conduction by comparing its epicardial depolarization times to measurements (15) for pacing at the right ventricular apex (RVA) and left ventricular (LV) free wall. In the present study we investigated whether the heterogeneous distribution of local amplitude and time course of contraction and reduced pump function during RVA pacing can be simulated realistically by the same model. To that purpose, simulated local circumferential strain in the LV midwall was compared with measured strain, obtained in a MRI tagging experiment. In particular, we focused on strain tracings in early, middle, and late depolarized regions and on the relation between midwall circumferential shortening during ejection and epicardial depolarization time. Results from the present study could be helpful in designing new ways of ventricular pacing.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
Animal protocols were approved by the National Heart, Lung, and Blood Institute Animal Care and Use Committee.
Maps of subepicardial depolarization timing and circumferential midwall strain for RVA pacing were measured in three dogs as described previously (8). In brief, a sock with 128 electrodes was placed over the ventricular epicardium of an anesthetized dog. Anesthesia was induced with an initial intravenous injection of thiopental sodium (25 mg/ml at 0.5 ml/kg) and maintained after endotracheal intubation with isoflurane (0.8–2%; Siemens ventilator 900D). Bipolar epicardial pacing electrodes were sewn onto the right atrium in all dogs. Similar electrodes were sewn onto the right ventricular (RV) apical free wall. RV pacing was done at a pacing rate of 110–125 beats/min, 
10–20% above the intrinsic rate. Pacing current was set to 20% above that needed for capture. Intrinsic sinus activation was suppressed by simultaneously pacing the right atrium. The animals were positioned in a MR scanner to obtain tagged cine images of the short and long axes of the heart during pacing. To view the onset of systole, the scanner trigger was delivered 40–50 ms before the pacing signal, initiating imaging in late diastole. The first image was captured 12 ms before the pacing signal.
Between image acquisitions, electrical recordings were obtained at an acquisition rate of 1,000 Hz. The animals were then euthanized, and the heart was excised. The coronary arteries were perfused from the aorta with isotonic saline at 50–60 mmHg to induce tissue turgor, and the heart was submerged in an isotonic saline bath to reduce body force deformation. Next, the saline in the LV and RV was displaced with vinyl polysiloxane (3M Express) injected through the corresponding atria and atrioventricular valves to fix the shape. The sock electrode locations were then recorded with a three-dimensional (3D) digitizer. From the tagged images, midwall circumferential strain was calculated throughout the cardiac cycle (20).
Unipolar voltage readings from each electrode were averaged over 20 heartbeats. Depolarization time tdep,epi was determined as the steepest downstroke of the electrode voltage reading. Depolarization time was related to the timing of the pacing trigger. Electrode locations in digitizer coordinates were transformed to scanner coordinates with a rigid-body rotation and translation.
Simulations
In the following paragraphs we give a short summary of the most important characteristics of the model of LV electromechanics, which have been reported previously (13, 15).
Geometry. The unloaded state was defined by zero transmural pressure. The LV wall was represented by a thick-walled, truncated prolate ellipsoid with a realistic myofiber orientation (Fig. 1A).
|
Hemodynamics. LV mitral inflow was simulated by increase of pressure in the nonactivated LV from 0 to 1 kPa. The LV was then activated by the pacing stimulus. LV pressure in the isovolumic contraction and relaxation phases was estimated such that LV cavity volume remained constant. The ejection phase started when LV cavity pressure exceeded aortic pressure, being set at 10 kPa. During ejection, LV cavity pressure was determined from the interaction of the LV with an aortic input impedance (40). Reversal of aortic flow ended the ejection phase.
Wall mechanics. Wall mechanics were determined by solving the equations of force equilibrium. Passive myocardial material was considered anisotropic and nonlinearly elastic (Fig. 1C). Active stress development in the myofibers depended on time, sarcomere length, and sarcomere shortening velocity (Fig. 1, D and E) and was initiated at the timing of depolarization.
Numerical implementation.
The eikonal-diffusion equation was solved by using a Galerkin finite-element method with eight-noded hexahedral elements with trilinear interpolation of the field of depolarization times. The LV wall was subdivided into 44,064 elements, resulting in 46,987 degrees of freedom and a spatial resolution of 1.4 mm.
The equations related to mechanics were solved by using a Galerkin finite-element method with 27-noded hexahedral elements with triquadratic interpolation of the displacement field. The LV wall was subdivided into 108 elements, with 3,213 degrees of freedom. All equations were solved on a 64-bit Origin 200 computer (SGI, Mountain View, CA), using a single processor at 225 MHz. The finite-element calculations were performed with the FORTRAN77 compiler-based package SEPRAN (SEPRA, Leidschendam, The Netherlands) on a UNIX platform. Calculation time for simulation of a complete cardiac cycle was 5 h.
Simulations and Data Analysis
Two separate simulations were performed. In the first simulation, normal sinus rhythm (SR) was simulated as in Ref. 13. Briefly, depolarization was started in three regions at the LV subendocardium below the equator and at one region at the septum at the RV side. At the subendocardium, prescribed wave velocity was six times higher than in the rest of the heart to account for a fully activated Purkinje system. Contraction was initiated with a heterogeneously distributed electromechanical delay such that mechanical activation was synchronous (13). Next, electromechanics after pacing at the RVA was simulated as in Ref. 15. Depolarization was started at the RVA. At the subendocardium, wave velocity was 1.7 times higher than in the rest of the heart. Contraction was initiated 0 ms after depolarization (Fig. 1B). The only difference in input for simulating mechanics on SR and RVA stimulations was the sequence of mechanical activation, while all other parameters were kept equal (geometry, passive and active material properties, preload, aortic impedance). For all simulations, global hemodynamics were computed as a function of time. In both experiment and simulation, circumferential natural strain 
cc in the midwall was defined as:
 | (1) |
Circumferential strain during ejection cc,ej was defined as the strain cc,eej at the end of ejection minus the strain cc,bej at the beginning of ejection:
 | (2) |
Algorithm to find beginning and end of ejection. To find the beginning of ejection in the experiment, we assumed that ventricular volume is closely related to mean of midwall circumferential strain, averaged over all regions of the LV wall, from which mean circumferential strain rate was computed (Fig. 2). From the largest shortening velocity (point A in Fig. 2), timing of 50% of the largest shortening velocity (point B) was determined. At this moment, the tangent line of shortening velocity (thin oblique line in Fig. 2) was computed. The intersection of the tangent line with zero shortening velocity yielded the beginning of ejection (point C). All midwall circumferential strains were set to 0 at the beginning of ejection, defining the reference state. The timing of minimum mean strain was chosen as the end of ejection (point D).
|
Early, middle, and late depolarized regions. We compared strain courses in more detail in early, middle, and late depolarized regions in the simulation and experiment for dog 3. Early, middle, and late depolarization were defined as 22–28%, 55–59%, and 94–100% of latest depolarization, respectively.
Relating ejection strain to depolarization time. It has been shown that local timing of depolarization determines myofiber strain during ejection (6), i.e., the local contribution to the ejected blood volume. We investigated this relation by computing the slope from a linear regression analysis (plus a 95% confidence interval) of the relation between epicardial depolarization time and midwall circumferential ejection strain in both experiment and simulation. The regression analysis was performed with the Analysis Toolpak in Microsoft Excel 2002. The three slopes from the experiment were averaged, yielding the mean ± SD slope. Thus midwall circumferential strain was assumed to closely match midwall myofiber strain, because myofibers are oriented within 20° from the circumferential direction (10). Also, epicardial depolarization time tdep,epi closest to the sites of calculated midwall circumferential strain was used to reproduce the methods in the experiment as closely as possible.
Stroke work density. In the experiment, circumferential strain was calculated from the midwall and depolarization at the epicardium. In the simulation, local depolarization times tdep (see Fig. 1B) and myofiber strain were calculated throughout the whole LV wall. Therefore, more information (and in addition stress and stroke work density) can be retrieved from the simulation than from the experiment performed.
Stroke work density Wf (J/m) was computed, using
 | (3) |
f and f represent total Cauchy myofiber stress and myofiber strain, respectively. |
RESULTS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
As shown in Fig. 3, maximum pressure and ejection fraction in the RVA pacing simulation (16.5 kPa and 48%, respectively) were lower than in the SR simulation (17.7 kPa and 55%, respectively). In the RVA simulation, ejection started later and peak aortic flow was lower than in the SR simulation.
|
In the SR simulation (Fig. 4, left), circumferential strain during ejection cc,ej was similar for all regions (–0.10 to –0.14). Both in the RVA simulation and measurements, strain patterns depended on epicardial depolarization time, as shown in Fig. 4, for early, middle, and late depolarized regions.
|
cc decreased during early isovolumic contraction. During ejection, cc increased initially. Thereafter, strain courses diverged such that total strain during ejection, cc,ej, ranged from slightly negative to positive (–0.02 to +0.05). Stretching continued in the isovolumic relaxation phase.
In the experiment the strain also decreased during early isovolumic contraction, reaching a minimum of about –0.05 30 ms before beginning of ejection. During the ejection phase, cc became positive up to 0.08 before midejection. Next, strain decreased slowly. Over the complete ejection phase, cc,ej was positive (0.00–0.06).
Middle depolarized regions.
In the RVA pacing simulation cc decreased by 0.08 in the isovolumic contraction phase. The decrease continued during ejection. Over the complete ejection phase, cc,ej was negative (–0.02 to –0.12).
In the experiment cc also changed by about –0.09 just before ejection. At end of ejection strain was negative (0 to –0.05).
Late depolarized regions.
In the RVA pacing simulation cc increased by 0.08 during the isovolumic contraction phase. Maximum cc occurred early in the ejection phase. Thereafter, strain decreased to values of –0.25 to –0.28 at end of ejection. Shortening continued in the beginning of isovolumic relaxation.
In the experiment cc increased by 0.10–0.12 during early isovolumic contraction, reaching a maximum of 0.03 30 ms before beginning of ejection. During ejection, cc,ej decreased to –0.12 to –0.27 at end of ejection. In the beginning of isovolumic relaxation, strain decrease continued.
Stress-Sarcomere Length Relation
Figure 5 shows myofiber stress as a function of sarcomere length (SL) in the RVA simulation in early, middle, and late depolarized regions. For the SR simulation myofiber stress and SL were taken from the same regions as in the RVA simulation. In the SR simulation the stress-SL loops progressed in a normal counterclockwise direction and were also similar in shape. Myofibers were stretched during filling (point A in Fig. 5). During isovolumic contraction (point B) some shortening was seen while stress was generated significantly. Stress reached a peak during ejection (point C), and myofibers shortened further until the isovolumic relaxation phase was reached (point D). During this last phase SL change was small and stress decreased until about zero at the end of the cycle.
|
Early depolarized region.
The stress-SL loop progressed in a clockwise direction. After filling (point A) myofibers shortened by 20% during isovolumic contraction (point B), while stress remained low, on the order of 1 kPa, in the first part of isovolumic contraction. In the second part of the isovolumic contraction phase, stress started to increase while the myofibers were stretched. Stretching continued during ejection (point C) and stress increased for the largest part of the ejection phase. During the first part of isovolumic relaxation (point D), SL still increased, until it decreased in the last part of this phase. Stress decreased throughout the whole isovolumic relaxation phase.
Middle depolarized region. The stress-SL loop progressed in a counterclockwise direction, as in the SR simulation. However, right after filling, myofibers were stretched first before they started to shorten while stress increased during isovolumic contraction. Both shortening during ejection and peak stress were smaller.
Late depolarized region. Also in this region, the stress-SL loop progressed in a counterclockwise direction, but here both SL and stress increased throughout the isovolumic contraction phase. Shortening during ejection and peak stress were both larger than in the SR simulation. Shortening also continued during the first part of isovolumic relaxation.
The area of myofiber stress-strain (replacing sarcomere length with strain in Fig. 5) loops represents stroke work density Wf (kJ/m), as computed from Eq. 3. Figure 6 shows a 3D representation of the LV wall, in which stroke work density for both simulations is presented in color code. Comparison with Fig. 1B shows that stroke work density was dependent on depolarization time as well. In the RVA simulation, it gradually increased from negative values (less than –3.0 kJ/m3) in early depolarized to supranormal values in late depolarized regions (>10 kJ/m3). Mean ± SD work in the RVA and SR simulation was 3.74 ± 4.0 and 4.90 ± 0.8 kJ/m3, respectively.
|
cc,ej is plotted as a function of epicardial depolarization time tdep,epi for the RVA simulation and the three dogs.
|
cc,ej and epicardial depolarization time tdep,epi for the RVA pacing simulation (–3.80 ± 0.28 s–1), was similar to the mean slope for the experiments [–3.34 s–1 (SD 0.97); n = 3]. |
DISCUSSION |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
Model
In the model geometry has been simplified. The real cardiac geometry with a LV and a RV is more complex than the truncated ellipsoid we used in the simulations. The load the RV exerts on the LV wall through RV pressure and the force transmission in the LV-RV attachment were neglected. During RVA pacing, RV pressure rises earlier than LV pressure, as observed in patients with left bundle branch block (11). The effect of these geometric simplifications might be investigated in a future study.
The inotropic states in the SR and RVA simulations were identical. Because this may vary from animal to animal, the effects of inotropy also might be investigated in more detail in a future study.
Choice of Reference State
To compare strains in simulation and experiments, an identical reference state must be chosen. Our choice was determined by considerations with respect to the experiment. Often, timing of the beginning or end of the ejection phase is determined by an aortic flow probe, but this approach is very difficult in the MRI scanner setup used. Alternatively, the LV pressure signal might be used. We considered this not very appropriate, because determining end diastole or beginning of ejection from LV pressure signals may be subject to larger error due to the gradual onset of LV pressure rise (end diastole) and the absence of a clear transition (beginning of ejection). We have also attempted to use the timing of stimulus as a reference state. However, in the heart there is a delay from excitation to contraction, so that time of electrical stimulation precedes true end diastole by a few tens of milliseconds. Therefore, we assumed that ventricular volume is closely related to the mean value of circumferential strain in the entire LV midwall and derived the moment of beginning ejection from the mean strain rate signal.
Measurements of Depolarization
Depolarization was measured at the epicardium. Previously performed studies have mapped endocardial depolarization (7, 34, 36). However, epicardial maps show a wider range of depolarization times than endocardial maps, which makes the calculation of the slope of the linear fit between myofiber ejection strain and depolarization time less prone to error. Ideally, one should measure depolarization and strain at the same location, but this is not feasible because of the limited accuracy of strain measurement at endocardium and epicardium and the requirement of invasive techniques to measure electrical activation at multiple sites in the midmyocardium.
Results from Simulations and Experiment
In the SR simulation, epicardial breakthrough occurred near the apex. After breakthrough, epicardial depolarization spread smoothly from apex to base, contrary to the finding of Arisi et al. (3). In the latter study, many epicardial breakthrough regions were measured in a mosaiclike pattern. This group addressed the occurrence of many breakthrough regions to bulges of excitation wave fronts that emerge from the trabecular endocardial surface and to discontinuous junctions between the Purkinje system and myocardium. In our model, the endocardium was smooth and we assumed a continuous layer of fast propagation at the endocardium. However, most characteristics (early activated regions, direction of propagation, timing) of the simulated normal propagating wave are physiological, as we discussed previously (13).
Time courses of circumferential strain cc are similar in the experiment and the RVA simulation (Fig. 4). The findings from the experiment were also similar to previously published experimental studies (4, 38, 41).
The main difference is related to timing, despite an identical procedure to determine the moment of beginning ejection in model and experiment. In early depolarized regions, minimum strain is reached earlier in the experiment (20 ms before beginning of ejection) than in the model (at beginning of ejection). Similarly, in late depolarized regions, maximum strain is reached earlier in the experiment (20 ms before beginning of ejection) than in the model (10 ms after beginning of ejection). To investigate the influence of the definition of beginning of ejection on the slope of the linear fit between myofiber ejection strain and depolarization time, we shifted the beginning of ejection in the simulation. Shifting the beginning of ejection such that timing of minimum strain in the early depolarized region in the simulation is equal to the timing in the experiment results in a slope of –3.19 ± 0.22 s–1. Shifting such that timing of maximum strain in the late depolarized region in the simulation is equal to the timing in the experiment results in a slope of –2.78 ± 0.18 s–1. These slopes are still in the range of the experiment.
Experimentally it has been made plausible that the relation between activation time and systolic shortening reflects the enhancement of late systolic contraction by early systolic stretch. The fact that the model was able to estimate the slope of this relation quite accurately indicates that the complex myocardial strain patterns occurring during ventricular pacing can be regarded as being based on differences in early and late depolarized regions. From the RVA simulation it was observed that in early depolarized regions myofibers shorten almost unloaded (Fig. 5). That is, cavity pressure is still low and the rest of the ventricular wall is not activated yet. Because all valves are closed, myofibers in the rest of the wall are stretched when early depolarized myofibers contract. When more myofibers become depolarized, stress in the wall increases and cavity pressure starts to rise. Early depolarized myofibers also become more loaded and start to generate stress.
Myofibers in late depolarized regions are stretched up to 10% and hence generate more stress than early depolarized myofibers because of the SL-active stress relation (Fig. 1D), the basis for the Frank-Starling mechanism. Therefore, during ejection, these late depolarized myofibers are able to shorten by >25% while stretching the early depolarized myofibers. Latter myofibers are stretched because of their lower stress, because their strain is low and they have reached the end of their twitch (they were activated 100 ms before myofibers in late depolarized regions). Hence, late depolarized myofibers contribute more to the ejection of blood than the rest of the heart (and also more than myofibers in a normal heart, as seen from the SR simulation), whereas the early depolarized ones even counteract that by being stretched.
The pattern of stroke work density Wf in the RVA simulation is similar to that reported earlier in the literature (23), with negative work in early depolarized regions and large positive work in late depolarized regions. Also, stroke work density is closely related to regional myocardial blood flow. Lindner et al. (17) showed that in patients with left bundle branch block myocardial blood flow is increased in late depolarized regions. The increase of Wf with depolarization time is explained by the fact that in early depolarized regions myofibers shorten in the isovolumic phase and are stretched during ejection and isovolumic relaxation. In late depolarized regions, SL at beginning of ejection is high, which enables large shortening and stress during ejection (and hence in need of more oxygen and blood flow). Overall, mean Wf is lower during ventricular pacing, compared with a normal heartbeat, which is also reflected in a smaller pressure-volume loop (Fig. 3).
In the present simulation of mechanics for normal SR a fixed, albeit regionally different, electromechanical delay was assumed, which fits with the data from the experiments in nonfailing dog hearts. It is being questioned whether electromechanical coupling in normal and failing hearts is the same. Leclercq et al. (16) found that in dyssynchronous failing hearts LV pacing improved mechanical synchrony and function without reducing electrical dyssynchrony. It is not clear from their data whether this result is due to altered local electromechanical delay or other factors like diastolic ventricular coupling (5). To simulate the experiment of Leclercq et al., it therefore seems necessary to use a biventricular model (14, 30), to account for ventricular interaction. Also, incorporation of the many changes at the cellular, tissue, and organ levels that are known to occur in failing hearts (29) will be required.
Usyk and McCulloch (32) concluded that the regional sequence of onset of myofiber shortening is an unreliable surrogate for regional depolarization. In our study, however, the goal was not to find a reliable surrogate for regional depolarization. We have focused on myofiber strain over the ejection phase in relation to depolarization and its agreement with experiments.
The resulting asynchrony of myofiber contraction during ventricular pacing affects pump function (21). In the long term, myocardial tissue structure changes (33) and may even contribute to the development of heart failure (2). Therefore, investigators are searching for better sites of pacing for optimal pump function. Positioning of the pacing electrode by trial and error is cumbersome. In assessing various pacing sites for minimal mechanical asynchrony, our mathematical model of cardiac electromechanics in the ventricle is likely to be a useful tool.
In conclusion, applying known principles of propagation of depolarization, time-dependent contraction of myofibers, and equilibria of forces, a finite-element model of LV electromechanics was developed. With this model we could realistically simulate mechanics in a ventricle that is paced at the RVA by only changing the sequence of depolarization, showing that the apparently complex events observed during pacing are caused by well-known basic physiological processes.
|
GRANTS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
|
FOOTNOTES |
|---|
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
|
REFERENCES |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSREFERENCES |
|---|
This article has been cited by other articles:
|
|
 |
|
  B. A. Coppola, J. W. Covell, A. D. McCulloch, and J. H. Omens Asynchrony of ventricular activation affects magnitude and timing of fiber stretch in late-activated regions of the canine heart Am J Physiol Heart Circ Physiol, July 1, 2007; 293(1): H754 - H761. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
N. H. L. Kuijpers, H. M. M. ten Eikelder, P. H. M. Bovendeerd, S. Verheule, T. Arts, and P. A. J. Hilbers Mechanoelectric feedback leads to conduction slowing and block in acutely dilated atria: a modeling study of cardiac electromechanics Am J Physiol Heart Circ Physiol, June 1, 2007; 292(6): H2832 - H2853. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
and 
, Opening and closing of valves for the SR and RVA simulations, respectively. Arrows indicate direction of the loop progression through time. Note that in the early depolarized region the loop for the RVA simulation progresses in a clockwise direction, whereas all other loops progress in a counterclockwise direction.
