Home | About | Journals | Submit | Contact Us | Français |

**|**Springer Open Choice**|**PMC2940053

Formats

Article sections

Authors

Related links

Annals of Biomedical Engineering

Ann Biomed Eng. 2010 October; 38(10): 3112–3123.

Published online 2010 May 25. doi: 10.1007/s10439-010-0071-x

PMCID: PMC2940053

NIHMSID: NIHMS222985

Niels F. Otani,^{}^{1} Stefan Luther,^{1,}^{2,}^{3} Rupinder Singh,^{4} and Robert F. Gilmour, Jr.^{1}

Niels F. Otani, Email: ude.llenroc@1ofn.

Associate Editor Anne Clough oversaw the review of this article.

Received 2010 March 23; Accepted 2010 May 10.

Copyright © The Author(s) 2010

The pattern of action potential propagation during various tachyarrhythmias is strongly suspected to be composed of multiple re-entrant waves, but has never been imaged in detail deep within myocardial tissue. An understanding of the nature and dynamics of these waves is important in the development of appropriate electrical or pharmacological treatments for these pathological conditions. We propose a new imaging modality that uses ultrasound to visualize the patterns of propagation of these waves through the mechanical deformations they induce. The new method would have the distinct advantage of being able to visualize these waves deep within cardiac tissue. In this article, we describe one step that would be necessary in this imaging process—the conversion of these deformations into the action potential induced active stresses that produced them. We demonstrate that, because the active stress induced by an action potential is, to a good approximation, only nonzero along the local fiber direction, the problem in our case is actually overdetermined, allowing us to obtain a complete solution. Use of two- rather than three-dimensional displacement data, noise in these displacements, and/or errors in the measurements of the fiber orientations all produce substantial but acceptable errors in the solution. We conclude that the reconstruction of action potential-induced active stress from the deformation it causes appears possible, and that, therefore, the path is open to the development of the new imaging modality.

Current technology does not allow us to measure large-scale action potential propagation deep inside the walls of the heart. This deficiency has a number of consequences. In the clinical setting, we rely on diagnostics such as the electrocardiogram (ECG) or electrophysiology (EP) studies to provide information about the heart’s electrical state. While these tests do certainly yield a wealth of valuable information, they often do not provide a fundamental understanding of the *dynamics* underlying the patient’s cardiac rhythm. For example, the field of cardiology is currently moving away from EP studies as a diagnostic for determining whether to implant an implantable cardioverter defibrillator (ICD),6 relying instead on the ejection fraction,1,9,20 paradoxically a mechanical rather than electrical measure, as the basis for the decision. This suggests that an improved diagnostic for assessing the patient’s rhythm dynamics is needed. At the research level, this lack of a deep, panoramic view of action potential activity means that we cannot even say with certainty that ventricular fibrillation (VF) is caused by multiple re-entrant action potential waves, although of course it is strongly suspected. Even if VF is a manifestation of multiple waves, the nature of the dynamics of these waves remains controversial, with their induction and maintenance being attributed to either tissue heterogeneity27,33 or steep electrical restitution,11,12,17 and their spatial patterning thought to be composed of either several rotating waves on equal dynamical footing,18,19,22,30 or one dominant “mother rotor” wave driving fibrillatory conduction.4,26 The lack of clarity on these issues has greatly complicated the development of effective therapies for the prevention and treatment of several types of tachyarrhythmias, including VF.

This bewildering landscape of theoretical wave propagation patterns and the mechanisms for their initiation and maintenance are what have led researchers to develop tools that might be used to actually determine what patterns exist during VF and other rapid, abnormal rhythms. Rudy *et al.*25,29 and others have developed a noninvasive method called ECGi to track action potential propagation on the epicardium using hundreds of electrodes situated on the body surface. This method does not as yet provide transmural information, however. Other approaches, such as plunge electrodes,8,10,24 transillumination,2 and optical tomography15,31 can provide a 3D representation of the electrical activity in the heart. The spatial resolution of the plunge electrode method is, however, relatively low. Furthermore, it is likely that the presence of these electrodes creates electrical heterogeneities that strongly modify the patterns of propagation.3 Transillumination and optical tomography are promising optical techniques; however, they require the use of voltage sensitive dyes that are phototoxic and prevent the application of these techniques *in vivo*. Furthermore, optical tomography cannot achieve sufficient temporal resolution to resolve transient phenomena present during fibrillation.

Magnetic resonance imaging (MRI) and computed tomography (CT) are also capable of noninvasively mapping the 3D contractile motion resulting from electrical activity in the heart. Although these imaging modalities provide very high spatial resolution, neither technology has the time resolution necessary (5–10 ms) to visualize, for non-repeating events, the relatively fast timescales inherent in the propagation of action potential waves.

We can argue that the 3D implementations of ultrasound are among the most promising for investigating the behavior of arrhythmias deep within cardiac tissue. Compared to the aforementioned imaging techniques, ultrasound imaging is attractive because of its relative ease of use, real-time feedback, portability, absence of ionizing radiation, cost effectiveness, and widespread availability in research and clinical applications. We can estimate the required temporal and spatial resolutions to be on the order of 10 ms and 1.5 mm, respectively. These are based on capturing 10 frames per rotation of a spiral wave rotating at 10 Hz, and the characteristic length scale associated with the calcium dynamics (a 10 ms timescale on a wave moving at 15 cm/s). Existing 2D ultrasound machines are easily capable of such resolutions. Indeed, ultrasound techniques have been demonstrated, at least in 2D, to be capable of yielding real-time, highly temporally resolved images of electromechanical wave contractions in humans.23,28

The proposed new imaging modality is based on the ability of ultrasound to detect wave-induced tissue deformation at depth within myocardial tissue. The idea is to extend conventional ultrasound strain imaging, which extracts strain information from the deformations appearing in the ultrasound images that are caused by action potential propagation activity, and convert these strains into the action potential-induced stresses that caused them. These stresses will then serve as markers for the presence of the action potentials, allowing us to create mappings of their propagation as functions of both time and space. A similar ultrasound technique called harmonic motion imaging is already the subject of intense study by the Konofagou and collaborators,16 although their goals, which are principally to assess the normal and pathological states of the cardiac tissue, are quite different from the action potential mapping goals of our research.

Here we report on a mathematically based computational technique that is capable of reconstructing the apparent motion of active stresses induced by the passage of action potentials from the mechanical tissue motion. We will regard this motion as motion that can be extracted from an ultrasound image sequence. The locations of these active stresses as functions of time will then, in many situations, serve as reasonable markers for the locations of the action potentials themselves, thereby allowing the study of propagation patterns of action potentials within cardiac tissue, in two and three spatial dimensions.

Our proposed imaging technique is based on the fact that propagating action potentials in the heart create calcium transients within the myocytes that induce active contractile stresses, as depicted by the left panel in Fig. 1. These active stresses are oriented predominantly in the local fiber direction (black arrows in the center panel), in the portion of the tissue occupied or recently occupied by the action potential. The active stresses, in turn, lead in general to the deformation of the entire tissue as illustrated schematically in the right panel of Fig. 1. Our method works by extracting the tissue displacements from ultrasound images showing the tissue deformation and converting them into the action potential-induced stresses that caused them, using a numerical algorithm we refer to here as the *inverse model* (as shown in Fig. 1). As part of the process we used to test the inverse model algorithm, we also built a *forward model* that simulates the biomechanical response of the 3D cardiac tissue to a given time-dependent active stress field meant to represent the stress created by a propagating action potential.

As indicated by the red arrows, the passage of action potentials through individual cells causes time-dependent processes in these cells (left panel), which lead to local contraction of the tissue primarily in the fiber direction (black arrows in the **...**

*Mathematical foundation for both the forward and inverse models.* Both models were based on the following standard equations:

1

and

2

Here *X*_{N} and *x*_{j} (*N*, *j* = 1, 2, 3), are the coordinates in the undeformed and deformed coordinate systems, respectively, *T*_{MN}(*E*) is the passive second Piola–Kirchhoff stress tensor, which depends on the Lagrange–Green strain tensor where the Cauchy–Green deformation tensor and the deformation gradient tensor *p* is the local hydrostatic pressure, and is the active second Piola–Kirchhoff stress tensor due to the action potential-induced force generated parallel to the orientation of the myocardial fibers. Sums over repeated indices are implicit. Equation (1) is a statement of force balance and ignores inertia, a standard assumption for the timescales we are considering.32 It also assumes the myocardium is a nonlinear, elastic medium. While the myocardium is perhaps better described as pseudoelastic, the variations of the stretch ratios (strains) for any given stress are no more than 20%,7 so we find this simplification appropriate. In this article, we use the hyperelastic, orthotropic pole-zero model of Nash and Hunter21 for the function *T*_{MN}(*E*). Equation (2) together with the appearance of the hydrostatic pressure in Eq. (1) expresses the assumption of incompressibility, common in the study of cardiac biomechanics.5,14,21

For our first attempt in studying this problem, we use the linearized version of Eqs. (1) and (2),

3

and

4

where, , *p* and the displacements are assumed small (i.e., *O*(*ε*), where *ε* is a small number). The indices *M*, *N*, *P*, and *Q* run over the values 1, 2, and 3, with summation again implied over repeated indices. This approximation should also be valid during ventricular fibrillation when the tissue deformations are typically very small; however, for cases of sinus rhythm or external electrical pacing, we expect we will have to use the fully nonlinear equations (1) and (2).

Finite element versions of Eqs. (3) and (4) were generated for the case of a tissue whose undeformed shape was that of a three-dimensional rectangular solid. Rectangular elements were defined in the obvious way relative to a 3D rectangular lattice of nodes whose edges coincided with the edges of the undeformed tissue. The spacings between nodal points, denoted Δ*X*, Δ*Y*, and Δ*Z*, were chosen to be uniform in each direction. These spacings may be arbitrarily chosen, as nothing in the resulting equations depends on their values. Weighting functions for Eqs. (3) took the form for a given node *l*, where the sum is over adjacent elements *e* and The element shape function is defined as when the node *l* is located at and the element *e* is with all other element shape functions defined in analogous fashion. The displacements were defined using trial functions identical to these weighting functions; i.e., where the ’s are constants. In contrast, since *p*, , and appear in Eq. (3) with one less derivative than the displacements, it is appropriate to define them in terms of amplitudes defined on the elements rather than the nodes; e.g., etc. Here we use a lower order trial function, defined simply as for **X** within the element *e*, and 0 otherwise. This lower order function was also used as the weighting function for Eq. (4). When these finite element representations were substituted into Eqs. (3) and (4), multiplied by the weighting functions, and integrated over **X**, we obtained,

5

6

Integration by parts was performed in deriving Eq. (5), assuming homogeneous Neumann boundary conditions (i.e., force-free conditions on all the surfaces). We also assume that the active second Piola–Kirchhoff stress tensor has only one nonzero component, *t*^{active}, in the 1–1 direction, before rotation by a unitary matrix *U* orients the tensor in the local fiber direction. This corresponds to our assumption that, as an action potential propagates by in some arbitrary direction, active contractile stress can only occur in the fiber direction, which we define to be in the “1” direction before rotation. This assumption is not trivial physiologically, as the orientation of fibers along a clear direction can, under some pathological conditions, fail to exist. Sums over *e* and *f* are over elements (only a single sum in the middle term), sums over *l* are over nodes, and sums over *M*, *P*, and *Q* are over 1, 2, and 3. Equation (5) depends on *N* = 1, 2, 3 and node index *j*; Eq. (6) just depends on *e*. Equations (5) and (6) are linear and homogeneous in the displacements defined at the nodes, the active stresses along the local fiber direction in each of the elements, , and the pressures in the elements, *p*^{e}. We note that any arbitrary rigid rotation or translation of a solution to Eqs. (5) and (6) is again a solution; thus it was useful to constrain the solution to a particular translation and rotational orientation by adding six equations to Eqs. (5) and (6): three that hold the mean displacement of the nodes to 0 in all three dimensions (which also translates into no net displacement of the center of mass if uniform mass density can be assumed):

7

and three equations that prevent rotation around the center of mass:

8

where The equations therefore have the form,

9

where all the elements in the first “super” matrix on the left are themselves matrices, and all the elements of the second matrix are column vectors. The first three rows of the matrix on the left of Eq. (9) correspond to the three components of Eq. (5), the fourth row corresponds to Eq. (6), and the fifth row imposes the constraints defined by Eqs. (7) and (8).

*The “forward” model: Determination of the tissue displacements caused by a given active stress field.* This function is easily accomplished by rewriting Eq. (9) as,

10

Given values for all the active stresses in the column vector *t*^{active}, we have four sets of equations (i.e., those associated with the first four rows of the matrix on the left) for the four unknown column vectors δ*x*_{1}, δ*x*_{2}, δ*x*_{3}, and *p*. These four sets of equations, when combined with the six constraints defined by the fifth row, have sufficient rank to allow solution using the built-in sparse matrix equation solver in the Matlab programming environment (The Mathworks, Inc.), thus yielding a complete description of the tissue deformation and compensating hydrostatic pressure produced by a given active stress field. Since there are no time derivatives in these equations, the temporal evolution of the deformations induced by a time-dependent active stress field may be created simply by solving Eq. (10) repeatedly for each instant in time for which we have active stress field data.

*The inverse model: Determination of the active stresses produced by action potentials using tissue displacement data.* Two variants of the inverse model were tested. Variant 1: If two of the three direction components (say, the “1” and “3” components) of the displacements can be extracted from the full three-dimensional volume, the model can reconstruct the active stress field that caused the displacements everywhere within the volume. We expect that this method will be useful in future experiments that make use of 3D ultrasound imaging techniques. The appropriate form of Eq. (9) for this method is,

11

In this case, the four sets of equations plus the six constraints represented by Eq. (11) actually constitute an overdetermined problem for the three unknown fields, δ*x*_{2}, *p*, and *t*^{active}. We therefore solved for them in the least-squares sense, again using Matlab’s built-in sparse matrix solvers.

Variant 2: If the deformation data come from a two-dimensional imaging plane, say, a plane oriented in the *x*_{1} and *x*_{3} (i.e., the *x* and *z*) directions, embedded within the volume, and consist of the two in-plane directional components of the displacement, then we can estimate, within the same plane, the active stress field that caused it. The solution to this problem is, in principle, relevant to the two-dimensional B-mode ultrasound imaging techniques. In this case, however, there is technically not enough data to uniquely determine the stress field. The problem is that the spatial derivatives in the direction perpendicular to the plane contribute to the determination of the stress. These derivatives cannot be determined solely from displacement data obtained from the plane. One of these derivatives may be obtained from the incompressibility condition; i.e., however, this still leaves the other two normal derivatives, and undetermined. For this initial study, we make the assumption that all planes parallel to this plane contain identical data for the displacement fields δ*x*_{1} and δ*x*_{3}. We then apply Eq. (11). We can expect this assumption to be valid if the variation of mechanical deformation is small in the *x*_{2} (i.e., the *y*) direction.

As with the forward problem, either inverse model can generate the active stress field by simply solving Eq. (11) for each time for which displacement data are available.

To test our methods, the forward model was used to generate the mechanical deformations induced by the traveling plane wave of active stresses shown in Fig. 2a in the cubic block of tissue shown Fig. 2b. The plane wave was chosen to have a width that was five-ninths the edge length of the cube, and contained a constant active stress of magnitude 2.0. This wave was used to mimic the contractile force that would be produced by an action potential plane wave. The block of tissue contained fibers whose direction rotated with depth, as is also the case in actual myocardial tissue. The fibers were oriented in the *y*-direction on the top surface of the cube-shaped tissue in Fig. 2b, rotated in the clockwise direction as a linear function of *z* to point in the *x*-direction in the horizontal plane in the middle of the cube, and continued to rotate linearly clockwise so that they pointed back in the *y*-direction at the bottom of the cube. For this study, a 10 × 10 × 10 lattice of nodes was employed.

(a) Applied active stresses in the local fiber direction as a function of *x* and *z* for all values of *y* at three selected times. The applied active stress was chosen to be in the gray regions, where is the local fiber direction, and zero everywhere else. **...**

We then used the deformations obtained from the forward model as input into the two variants of the “inverse” model, based on Eq. (11). Variant 1, which takes as input two of the three displacement components expressed as spatially 3D data, was able to reproduce the original active stresses essentially exactly, to within 1 part in 10^{12}, as is readily apparent by comparing Figs. 2c to to2a.2a. Variant 2 was able to do a reasonable job of reconstructing the stress field plane wave in the plane from which the deformation data were taken (Fig. 2d). When random errors were purposely introduced into the fiber direction, both variants of the inverse code were still reasonably successful in reconstructing the active stresses, as illustrated in Figs. 2e and and2f,2f, suggesting that the reconstruction process is not highly sensitive to errors we might make in measuring the fiber direction in the experiment. Here the errors were chosen from a distribution of zero mean and standard deviation of 10°.

Studies of this type were conducted for 10 different orientations of this rotated fiber pattern. Specifically, the ability of Variants 1 and 2 to reconstruct the active stress field of a plane wave propagating in the *x*-direction, in both the presence and absence of fiber direction error, was studied for the cases in which the entire fiber orientation spatial pattern just described was rotated through angles ranging from 0° to 90° around the *z*-axis in increments of 10°. For all 10 orientations, Variant 1 without fiber error reproduced the original active stress field to within 1 part in 10^{12}. For the other three cases (i.e., Variant 2 without fiber errors and both variants with fiber errors), we found that the reconstruction was more accurate in regions in which the fiber orientations lie in the *x*–*z* plane. For example, for the fiber orientation pattern used in all the plots in Fig. 2, this region intersects the displayed cross section in the white rectangle in each panel of Figs. 2d–2f. The higher accuracy within each rectangle is manifested as more red, orange, and yellow cells in the portion of the rectangle inside original non-zero active stress region, and primarily blue cells outside this region (but still inside the rectangle). The position of the region containing *x*–*z* oriented fibers changes in each of the 10 orientations, but in each case, the accuracy was found to be highest within this region. Other than differences in the positioning of this higher accuracy region, the reconstructed active stresses for the 10 cases did not exhibit features substantively different from those shown in Figs. 2c–2f.

When the magnitudes of the reconstructed active stresses of all the cells inside (red) and outside (blue) the original non-zero active stress regions at all measured times from all 10 of these studies were plotted as distributions, we obtained the plots shown in Fig. 3. As expected, for Variant 1 with no fiber orientation errors, there is no spread of the values around the original values of 2.0 and 0.0, respectively. (The slight amount of spread shown in Fig. 3a is entirely due to the width of the bins used to create the distributions.) The distributions for the other three cases showed substantial spread, with somewhat more spread exhibited by the cells inside the original non-zero active stress region for the Variant 2 cases (Figs. 2b and and2d)2d) than the Variant 1 case (Fig. 2c).

Distribution of magnitudes of the reconstructed active stresses for all cells inside (red) and outside (blue) the nonzero active stress region, as calculated by (a) Variant 1 of the inverse model, (b) Variant 2, (c, d) Variants 1 and 2, respectively, **...**

Despite this spread, the two distributions are fairly well separated from each other. Thus, a value, , may be chosen such that, if any given cell has it is most likely in the active stress region, while if it is mostly likely not. Figure 4 shows how well this classification scheme works for different values of

Fraction of cells misclassified as either inside the non-zero active stress region when they are classified as being outside, or vice versa, as a function of the classification parameter . All cells from all recorded times from all 10 of the runs described **...**

We observe that, for the Variant 1 of the inverse model (i.e., the variant that uses 3D data as input), the fraction of cells misclassified is smallest when , while for Variant 2, the optimal value for is 0.63 or 0.61, depending on whether fiber error is present or not. Using as the classification criterion for the data shown in the two Variant 1 runs shown in Fig. 2 (i.e., panels c and e), we obtain Figs. 2g and and2i,2i, respectively, where orange represents cells in classified as being inside the active stress region, and blue outside. The agreement is perfect between Figs. 2g and and2a,2a, as expected, since the solution is overdetermined in this case and not compromised by fiber orientation errors. The agreement between Figs. 2i and and2a,2a, when fiber error is present, is not perfect, but is acceptable, in the sense that it provides good insight into how the active stress region (and thus, presumably, the corresponding action potential) is propagating. When was used to classify the Variant 2 runs (e.g., Figs. 2d and and2f),2f), favorable results were again obtained. As shown in Figs. 2h and and2j,2j, a reasonably clear representation of the left-to-right motion of the active stress region is again obtained, despite the presence of a few misclassified cells.

Simulations were also conducted to develop insight into the effect that error and/or noise in the ultrasound image displacement data will have on our inverse model. Since we currently have no clear indication of the type of noise that might be present in these data, we used a white noise power spectrum for the error, which presumably is a fairly pessimistic assumption. Accordingly, normally distributed random errors with mean 0 and standard deviation 0.024 were added to each of the displacements calculated by the forward model for the cubic system. This value for the standard deviation represented an error of 10% in the displacement measurements, since the displacements calculated from the forward model were found to have a spread in magnitude characterized by a standard deviation of 0.24. The effects of these errors for Variant 1 of the inverse model with no fiber errors are shown in Fig. 5.

As seen in Fig. 5a, there is significant spread in the calculated values of the active stress, both in the active stress “wave” (red curve) and outside it (blue curve). However, if we classify the pixels as being in the wave when the calculated active stress is greater than 1, and outside the wave otherwise, we found that 8.2% of all the pixels were classified as in the AP when they were not, while 4.8% were classified as not in the AP when they were. The remaining 86.9% of the pixels were classified correctly. It was also interesting to note that the classification of pixels was better near the boundaries (e.g., Fig. 5c) than in the center (e.g., Fig. 5b).

Examination of the deformation of the tissue as shown in Fig. 2b suggests that the displacements created by a region of nonzero active tension are not confined to that region, but instead tend to appear throughout the tissue. This observation is supported by the following calculation, which derives the displacement field created by a point active stress. Suppose that the *zz* component of the active stress *T*_{zz} is only nonzero at the origin of a three-dimensional domain of infinite extent, while all other components of the active stress are 0 everywhere. Transformation of Eqs. (3) and (4) into a cylindrical coordinate system yields the following three equations:

12

13

14

where *δR* and *δZ* are displacements in the *R* and *Z* directions, respectively, and the constitutive relationship between strain and active stress is assumed to be isotropic with Young’s modulus *T*. We find that, when the following forms for the dynamical variables and *T*_{ZZ} are substituted into Eqs. (12–14):

15

16

17

18

algebraic expressions independent *R* and *Z* are obtained among the coefficients *C*_{R}, *C*_{Z}, *C*_{P}, and *C*_{T} indicating that Eqs. (15–18) are eigenmodes of the system of equations (12–14). Here *K* and *L* are arbitrary constants, and *J*_{0} and *J*_{1} are 0th and 1st order Bessel functions of the first kind. Rearrangement of these algebraic relationships allows all of these coefficients to be expressed in terms of *C*_{T}:

19

20

21

To apply these eigenmode solutions to the case of a point active stress at the origin, we note that an active stress field that is only nonzero at the origin having an integrated strength of *t*_{0} may be written as a Fourier–Bessel integral13:

22

where

23

Since the system of equations is linear, we can use the principal of superposition to write the displacements and hydrostatic pressure resulting from this point active stress as Fourier–Bessel integrals, with coefficients given by Eqs. (19–21), where *C*_{T} is given by Eq. (23). The resulting integral expressions may be evaluated using integral tables.13 Transforming the expressions from cylindrical to spherical (*r*,θ,) coordinates, we obtain the following:

24

25

26

These Green’s function solutions show that the displacements and hydrostatic pressure fall off as 1/*r*^{2} and 1/*r*^{3}, respectively, in response to active stress located at a single point. The response to active stress thus behaves like a long-range force not unlike gravity or the electrostatic force, affecting the displacements and pressure far from the actual source. This eliminates the possibility that a technique can simply use the local presence of displacements as a marker for the presence of action potentials activity nearby. Thus, it should not be possible in general to determine the location of the source, i.e., the active stresses, without considering the entire problem, including boundary conditions, within the framework of an electromechanical model. We also note that, while the displacements scale inversely with the Young’s Modulus *T* of the tissue, the pressure, which enforces the incompressibility condition, is independent of this quantity. Also, as expected, all quantities scale linearly with the magnitude of the active stress *t*_{0}.

In this article, we describe key mathematical techniques for visualizing the location and motion of waves of action potential-induced active stress traveling both on the surface and deep within myocardial tissue. As an initial test of these ideas, we have developed the linear version of the algorithm and studied its ability to reconstruct simple plane waves of active stress propagating in a three-dimensional cubic system containing fiber directions that rotate with depth. We find that, when the deformation data consist of two of the three components of the displacement field defined throughout the 3D system, the method is very precise, allowing essentially exact reconstruction of the active stresses. This result is to be expected—the theory shows that, when the active stresses are all a function of a single scalar field (i.e., the active stress along the local fiber direction) and the orientation of these fibers is known throughout the tissue, then the problem is actually overdetermined, and thus will yield the exact solution, when it exists. This calculation is the one relevant to ultrasound machines capable of imaging three-dimensional volumes within the myocardium.

In contrast, when the two components of displacement are only known within a two-dimensional imaging plane, as is the case for B-mode ultrasound imaging, determination of the active stresses within this same plane is an underdetermined problem, requiring reasonable assumptions to be made. As a first attempt, we assumed slow variation of the two in-plane displacements in the direction perpendicular to the imaging plane. With this assumption, we find that the method is still able to do a reliable reconstruction of a plane wave consisting of active stresses along the local fiber direction. This type of assumption may not be valid when the tissue properties are strongly heterogeneous, or when the pattern of wave propagation contains more spatial structure, as would be the case, for example, during ventricular fibrillation. This issue will be a focus of further investigation in future studies.

Even when errors of order 10° were purposely introduced into the fiber directions, both variants of the inverse model performed reasonably well in reconstructing the active stress plane wave. This result suggests that the reconstruction process is not highly sensitive to errors in the measurements of fiber orientation, an important characteristic of the method, as a requirement of precise measurements of the fiber orientations would greatly complicate the use of the method as an imaging modality.

The method also performed reasonably well when white noise of relative amplitude 10% was added to the displacement data. We currently have no predictions about the spatial power spectrum of noise and error we will find in the ultrasound data; however, we believe that the assumption of white noise should be a fairly pessimistic one. It is not too hard to imagine the existence of short-distance positive correlation in the error that results from a tendency for points in the tissue that are in close proximity to move together. Thus, we can reasonably expect that the performance of our inverse algorithm will be even better in practice than is shown here.

One can question whether such an involved inverse calculation is necessary to extract the active stress from the tissue displacements they induce. However, as we demonstrate in the calculations embodied by Eqs. (12–26), the magnitudes of these displacements fall off only as fast as 1/*r*^{2} with distance, so we believe it is unlikely that active stresses in a particular location can be determined solely from the displacements present in that area.

The method described here does not yield detailed membrane potential profiles of the action potentials. Such profiles are likely to be difficult or impossible to extract from the output of our method, the active stresses, because a clear one-to-one relationship between these stresses and the level of the membrane potential does not exist. This might, for example, be the case in partially refractory tissue or near the centers of spiral wave rotation. Nevertheless, we expect that the spatial and temporal patterning of active stresses yielded by our method should be a reasonable indicator of the location, motion, and spatial patterning of action potentials and their propagation, most of the time over most of the tissue. This level of imaging is, in our opinion, still likely to be very useful in determining the essential dynamics of the action potentials, which we believe will be important in diagnosing the nature of their propagation during clinically important phenomena such as ventricular polymorphic tachycardia and fibrillation.

While we are currently targeting ultrasound imaging as the application for this method, we note that the applicability of the method is not necessarily restricted to the ultrasound modality. Any imaging method that can be translated into mechanical tissue deformations in at least a 2D imaging plane is a potential candidate for use with this method.

S. Luther acknowledges support from the Bundesministerium für Bildung und Forschung (FKZ 01EZ0905/6). R. Singh was supported in part by the National Science Foundation, Grant No. 0800793. We thank VisualSonics Inc. (Toronto, Canada) for useful discussions. We appreciate the reviewers’ helpful and useful comments in preparing this article.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

1. Bardy GH, Lee KL, Mark DB, Poole JE, Packer DL, Boineau R, Domanski M, Troutman C, Anderson J, Johnson G, McNulty SE, Clapp-Channing N, Davidson-Ray LD, Fraulo ES, Fishbein DP, Luceri RM, Ip JH. Sudden Cardiac Death in Heart Failure Trial (SCD-HeFT) Investigators. Amiodarone or an implantable cardioverter-defibrillator for congestive heart failure. N. Engl. J. Med. 2005;352:225–237. doi: 10.1056/NEJMoa043399. [PubMed] [Cross Ref]

2. Baxter WT, Mironov SF, Zaitsev AV, Jalife J, Pertsov AM. Visualizing excitation waves inside cardiac muscle using transillumination. Biophys. J. 2001;80:516–530. doi: 10.1016/S0006-3495(01)76034-1. [PubMed] [Cross Ref]

3. Beaudoin DL, Roth BJ. Effect of plunge electrodes in active cardiac tissue with curving fibers. Heart Rhythm. 2004;4:476–481. [PubMed]

4. Chen J, Mandapati R, Berenfeld O, Skanes AC, Jalife J. High-frequency periodic sources underlie ventricular fibrillation in the isolated rabbit heart. Circ. Res. 2000;86(1):86–93. [PubMed]

5. Costa KD, Holmes JW, McCulloch AD. Modelling cardiac mechanical properties in three dimensions. Philos. Trans. Roy. Soc. Lond. 2001;A35:1233–1250. doi: 10.1098/rsta.2001.0828. [Cross Ref]

6. Daubert JP, Zareba W, Hall WJ, Schuger C, Corsello A, Leon AR, Andrews ML, McNitt S, Huang DT, Moss AJ. Predictive value of ventricular arrhythmia inducibility for subsequent ventricular tachycardia or ventricular fibrillation in Multicenter Automatic Defibrillator Implantation Trial (MADIT) II patients. J. Am. Coll. Cardiol. 2006;47:98–107. doi: 10.1016/j.jacc.2005.08.049. [PubMed] [Cross Ref]

7. Demer LL, Yin FCP. Passive biaxial mechanical properties of isolated canine myocardium. J. Physiol. 1983;339:615–630. [PubMed]

8. El-Sherif N, Chinushi M, Caref EB, Restivo M. Electrophysiological mechanism of the characteristic electrocardiographic morphology of torsade de pointes tachyarrhythmias in long-QT syndrome: detailed analysis of ventricular tridimensional activation patterns. Circulation. 1997;96:4392–4399. [PubMed]

9. Epstein A. E., J. P. DiMarco, K. A. Ellenbogen, N. A. M. Estes III, R. A. Freedman, L. S. Gettes, A. M. Gillinov, G. Gregoratos, S. C. Hammill, D. L. Hayes, M. A. Hlatky, L. K. Newby, R. L. Page, M. H. Schoenfeld, M. J. Silka, L. W. Stevenson, and M. O. Sweeney. ACC/AHA/HRS 2008 Guidelines for device-based therapy of cardiac rhythm abnormalities: a report of the American College of Cardiology/American Heart Association Task Force on Practice Guidelines (Writing Committee to Revise the ACC/AHA/NASPE 2002 Guideline Update for Implantation of Cardiac Pacemakers and Antiarrhythmia Devices): developed in collaboration with the American Association for Thoracic Surgery and Society of Thoracic Surgeons. *Circulation* 117:e350–e408, 2008. [PubMed]

10. Frazier DW, Wolf PD, Wharton JM, Tang AS, Smith WM, Ideker RE. Stimulus-induced critical point. Mechanism for electrical initiation of reentry in normal canine myocardium. J. Clin. Invest. 1989;83:1039–1052. doi: 10.1172/JCI113945. [PMC free article] [PubMed] [Cross Ref]

11. Garfinkel A, Kim YH, Voroshilovsky O, Qu Z, Kil JR, Lee MH, Karagueuzian HS, Weiss JN, Chen PS. Preventing ventricular fibrillation by flattening cardiac restitution. Proc. Natl Acad. Sci. USA. 2000;97:6061–6066. doi: 10.1073/pnas.090492697. [PubMed] [Cross Ref]

12. Gilmour RF, Jr, Chialvo DR. Electrical restitution, critical mass, and the riddle of fibrillation. J. Cardiovasc. Electrophysiol. 1999;10:1087–1089. doi: 10.1111/j.1540-8167.1999.tb00281.x. [PubMed] [Cross Ref]

13. Gradshteyn, I. S., and I. M. Ryzhik. Table of Integrals, Series, and Products: Seventh Edition, edited by A. Jeffrey and D. Zwillinger. Burlington, MA: Academic Press, 2007, pp. 659–753.

14. Guccione JM, McCulloch AD, Waldman LK. Passive material properties of intact ventricular myocardium determined from a cylindrical model. ASME J. Biomech. Eng. 1991;113:42–55. doi: 10.1115/1.2894084. [PubMed] [Cross Ref]

15. Hillman EMC, Bernus O, Pease E, Bouchard MB, Pertsov A. Depth-resolved optical imaging of transmural electrical propagation in perfused heart. Opt. Exp. 2007;15(26):17827–17841. doi: 10.1364/OE.15.017827. [PMC free article] [PubMed] [Cross Ref]

16. Hynynen, K. H., and E. Konofagou. Harmonic Motion Imaging. Patent 6,984,209, 10 Jan 2006.

17. Karma A. Electrical alternans and spiral wave breakup in cardiac tissue. Chaos. 1994;4:461–472. doi: 10.1063/1.166024. [PubMed] [Cross Ref]

18. Lee MH, Qu ZL, et al. Patterns of wave break during ventricular fibrillation in isolated swine right ventricle. Am. J. Physiol. Heart Circ. Physiol. 2001;281(1):H253–H265. [PubMed]

19. Moe GK. On the multiple wavelet hypothesis of atrial fibrillation. Arch. Int. Pharmacodyn. Ther. 1962;140:183–188.

20. Moss AJ, Zareba W, Hall WJ, Klein H, Wilber DJ, Cannom DS, Daubert JP, Higgins SL, Brown MW, Andrews ML. Multicenter Automatic Defibrillator Implantation Trial II Investigators. Prophylactic implantation of a defibrillator in patients with myocardial infarction and reduced ejection fraction. N. Engl. J. Med. 2002;346:877–883. doi: 10.1056/NEJMoa013474. [PubMed] [Cross Ref]

21. Nash MP, Hunter PJ. Computational mechanics of the heart. J. Elasticity. 2000;61:113–141. doi: 10.1023/A:1011084330767. [Cross Ref]

22. Panfilov A. Spiral breakup as a model of ventricular fibrillation. Chaos. 1998;8:57–64. doi: 10.1063/1.166287. [PubMed] [Cross Ref]

23. Pernot M, Fujikura K, Fung-Kee-Fung SD, Konofagou EE. ECG-gated, mechanical and electromechanical wave imaging of cardiovascular tissues in vivo. Ultrasound Med. Biol. 2007;33(7):1075–1085. doi: 10.1016/j.ultrasmedbio.2007.02.003. [PubMed] [Cross Ref]

24. Pogwizd SM, Corr PB. Reentrant and nonreentrant mechanisms contribute to arrhythmogenesis during early myocardial ischemia: results using three-dimensional mapping. Circ. Res. 1987;61:352–371. [PubMed]

25. Ramanathan C, Ghanem RN, Jia P, Ryu K, Rudy Y. Noninvasive electrocardiographic imaging for cardiac electrophysiology and arrhythmia. Nat. Med. 2004;10(4):422–428. doi: 10.1038/nm1011. [PMC free article] [PubMed] [Cross Ref]

26. Samie FH, Berenfeld O, Anumonwo J, Mironov SF, Udassi S, Beaumont J, Taffet S, Pertsov AM, Jalife J. Rectification of the background potassium current: a determinant of rotor dynamics in ventricular fibrillation. Circ. Res. 2001;89(12):1216–1223. doi: 10.1161/hh2401.100818. [PubMed] [Cross Ref]

27. Sampson KJ, Henriquez C. Simulation and prediction of functional block in the presence of structural and ionic heterogeneity. Am. J. Physiol. 2001;281:H2597–H2603. [PubMed]

28. Wang S, Lee WN, Provost J, Luo J, Konofagou EE. A composite high-frame-rate system for clinical cardiovascular imaging. IEEE Trans. Ultrason. Ferroelectr. Freq. Control. 2008;55:2221–2233. doi: 10.1109/TUFFC.921. [PubMed] [Cross Ref]

29. Wang Y, Rudy Y. Applications of the method of fundamental solutions to potential-based inverse electrocardiography. Ann. Biomed. Eng. 2006;34(8):1272–1288. doi: 10.1007/s10439-006-9131-7. [PMC free article] [PubMed] [Cross Ref]

30. Weiss JN, Chen PS, et al. Ventricular fibrillation—How do we stop the waves from breaking? Circ. Res. 2000;87(12):1103–1107. [PubMed]

31. Wellner M, Bernus O, Mironov SF, Pertsov AM. Multiplicative optical tomography of cardiac electrical activity. Phys. Med. Biol. 2006;51:4429–4446. doi: 10.1088/0031-9155/51/18/001. [PubMed] [Cross Ref]

32. Whiteley JP, Bishop MJ, Gavaghan DJ. Soft tissue modeling of cardiac fibres for use in coupled mechano-electric simulations. Bull. Math. Biol. 2007;69:2199–2225. doi: 10.1007/s11538-007-9213-1. [PubMed] [Cross Ref]

33. Yan GX, Shimizu W, Antzelevitch C. Characteristics and distribution of M cells in arterially perfused canine left ventricular wedge preparations. Circulation. 1998;98:1921–1927. [PubMed]

Articles from Springer Open Choice are provided here courtesy of **Springer**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |