A New Methodology for Multiscale Myocardial Deformation and Strain Analysis Based on Tagging MRI

Abstract

Myocardial deformation and strain can be investigated using suitably encoded cine MRI that admits disambiguation of material motion. Practical limitations currently restrict the analysis to in-plane motion in cross-sections of the heart (2D + time), but the proposed method readily generalizes to 3D + time. We propose a new, promising methodology, which departs from a multiscale algorithm that exploits local scale selection so as to obtain a robust estimate for the velocity gradient tensor field. Time evolution of the deformation tensor is governed by a first-order ordinary differential equation, which is completely determined by this velocity gradient tensor field. We solve this matrix-ODE analytically and present results obtained from healthy volunteers as well as from patient data. The proposed method requires only off-the-shelf algorithms and is readily applicable to planar or volumetric tagging MRI sampled on arbitrary coordinate grids.

1. Introduction

Cine MRI, combined with (C)SPAMM ((C)SPAMM = (Complementary) SPAtial Modulation of Magnetization) encoding technology [1–4], admits disambiguation of local tissue motion, thus enabling the extraction of myocardial deformation and strain [5], which are known to correlate with cardiac pathologies. In particular, Götte et al. found that strain is more accurate than geometry in discriminating dysfunctional from functional myocardium [6]. Deformation and strain can be operationalized in various ways, either without explicit a priori regularization, or through exploitation of sparse constraints combined with interpolation and/or regularization [5, 7–33]. Possible encodings are (DENSE) Displacement ENcoding with Stimulated Echoes, cf. Aletras et al. [2), and (HARP) HARmonic Phase, cf. Osman et al. [20, 21]. Tagging-based methods using HARP technology form our point of departure, compare with Figure 1 for an illustration.

(a)

(b)

Figure 1 

Vertical (a) and horizontal (b) SPAMM encoding in a basal short-axis slice at approximately mid-systole for a patient.

Given a dense motion field within the myocardium, our aim is to devise an operational procedure for direct extraction of myocardial deformation and strain. By “direct” we mean that we seek to obviate sophisticated preprocessing steps, such as segmentation of, or interpolation between tag lines, and finite element methods explicitly coupled to the tagging pattern. Although such sophisticated “indirect” procedures exist and have been proven powerful, they require specific algorithmics that is neither trivially implemented nor readily available. Instead we aim for a multiscale, optimally conditioned, intrinsically parallelizable, linear algorithm for obtaining deformation (and thus strain) in analytically closed form. Optimal conditioning is achieved by exploitation of the scale degree of freedom in the definition of spatiotemporal differential image structure. In addition we aim to minimize the number of extrinsic control parameters. We believe that the parsimony of our method facilitates applicability and optimization, since only off-the-shelf algorithms (linear filtering and inversion of linear systems) are needed in our computation of myocardial deformation and strain.

Our approach is as follows. To begin with, the velocity gradient tensor field, which is the prerequisite for the model we present in this paper, is obtained using the multiscale motion extraction algorithm for scalar image sequences proposed by Florack et al. [12]. This algorithm is adapted to the situation of tagging MRI data [10, 13]. Slick encoding and/or linear filtering yields an 𝑛-tuple of independent scalar phase images of some fixed spatiotemporal region of interest, in which 𝑛 is the dimension of space (𝑛=2,3, as may be the case). The intrinsic, 𝑛-fold underdetermined “optic flow constraint equation” can, by construction, be applied to each of these 𝑛 phase images, yielding an unambiguous system of equations for the underlying dense motion field to any desired differential order, and obviating Tikhonov regularization (which would bring in at least one additional control parameter) as a disambiguating prior (as is typically the case in optic flow applications due to the aperture problem). It has been verified that the first-order scheme indeed produces a plausible, dense, and robust velocity gradient tensor field within the myocardium. Details of its construction as well as the underlying automatic scale selection mechanism can be found in the cited literature (for an example of automatic scale selection, cf. also Niessen et al. [34]). In view of the cited work we will de-emphasize multiscale motion extraction and concentrate on subsequent deformation and strain analysis.

In Section 2 we outline in detail how to arrive at a closed-form analytical solution for the deformation tensor field (Section 2.1) and hence the strain tensor field (Section 2.2), given the velocity gradient tensor field. The novelty of this approach is that we circumvent numerical approximations in all intermediate steps, and introduce numerics only at the ultimate stage where we sample the resulting analytical tensor field expressions. Besides avoiding in this way numerical errors that may be difficult to quantify, this procedure has the advantage of being mathematically transparent, computationally trivial, and intrinsically parallellizable. (These properties, in fact, also hold for the multiscale motion extraction algorithm used to provide the input velocity gradient tensor field, loc. cit.) Some experimental results are given in Section 3 to demonstrate the feasibility of our approach.

Section 4 concludes our work with a summary and briefly sketches work in progress.

Data acquisition details for the scans used in this paper are given in the appendix.

2. Theory

2.1. Deformation

The velocity gradient tensor, with components 𝐿^𝛼_𝛽 relative to a coordinate frame, relates the rate of change of a momentary infinitesimal line element 𝑑𝑥̇^𝛼 to the line element 𝑑𝑥^𝛽 itself. (One can either view 𝑑𝑥^𝛼 as the 𝛼-component of an “infinitesimal vector” attached to a fiducial material point in Euclidean space and representing a directed material line element, that is, an “infinitesimal” number—or as the 𝛼th element of a local covector frame, depending on engineering or mathematical mind set.) From 𝑑𝑥̇^𝛼=𝑑𝑣^𝛼 it follows, using the chain rule, that

𝑑𝑥̇^𝛼=𝐿^𝛼_𝛽𝑑𝑥^𝛽with𝐿^𝛼_𝛽=𝜕𝑣^𝛼𝜕𝑥^𝛽(𝛼,𝛽=1,…,𝑛). (1) (The Einstein summation convention applied here will be used henceforth.) The numbers 𝐿^𝛼_𝛽 constitute a matrix 𝐋 with row and column index 𝛼,𝛽. To get an estimate of this tensor field we have applied the algorithm of Florack et al. and van Assen et al., solving a linear system of algebraic equations enforcing optimal conditioning through scale selection. The reader is referred to the literature for details [10, 12, 13].

Tissue deformation can be described by a map 𝑓∶𝑀→𝑁∶𝑥(𝑋,𝑡₀)↦𝑥(𝑋,𝑡), in which 𝑥(𝑋,𝑡) denotes the spatial position of a material point at time 𝑡, with reference position 𝑋=𝑥(𝑋,𝑡₀) at time 𝑡₀ (Lagrange picture). Domain 𝑀 and codomain 𝑁 are copies of the deformable tissue medium (a subset of ℝ^𝑛) at times 𝑡₀, respectively, 𝑡≥𝑡₀. We fix 𝑡₀ so as to correspond to end-diastole. (During the diastolic phase the heart relaxes and is refilled with blood. It is followed by the systolic phase, in which the heart contracts. It is this phase that traditionally receives most attention by cardiologists.)

The associated differential map,

𝑓_∗≡𝑑𝑓∶𝑇𝑀_𝑋→𝑇𝑁_𝑥∶𝑣↦𝑓_∗(𝑣)=𝑣^𝑖𝐹^𝛼_𝑖𝐟_𝛼, (2) relative to a basis {𝐟_𝛼} of 𝑇𝑁_𝑥, provides a local linearization of tissue deformation, that is, the deformation tensor. Here 𝑇𝑀_𝑋 is the tangent space of 𝑀 at the fiducial point 𝑋=𝑥(𝑋,𝑡₀), and 𝑇𝑁_𝑥 that of 𝑁 at the 𝑓-mapped point 𝑥=𝑥(𝑋,𝑡). The matrix 𝐅(𝑡,𝑡₀) of this linear map relative to the local coordinate charts {𝑥^𝛼,𝑁} and {𝑋^𝑖,𝑀} is given by the Jacobian

𝐹^𝛼_𝑖=𝜕𝑥^𝛼𝜕𝑋^𝑖. (3) By virtue of the chain rule, the relation between deformation and velocity gradient tensors, (1) and (3) is given by the first-order ODE [35]

𝐅̇=𝐋𝐅, (4) subject to an initial condition, namely, 𝐅(𝑡=𝑡₀,𝑡₀)=𝐈. Equivalently,

𝐅(𝑡,𝑡₀)=𝐈+^𝑡_𝑡0𝐋(𝑠)𝐅(𝑠,𝑡₀)𝑑𝑠. (5) Only for stationary flows, that is, 𝐋(𝑡)=𝐋₀ pointwise constant, this initial value problem admits a trivial solution 𝐅(𝑡,𝑡₀)=exp((𝑡−𝑡₀)𝐋₀). (For the sake of simplicity the spatial dependency of all field entities is suppressed in the notation.) However, we do not have stationarity, and so we must proceed more carefully.

Equation (5) induces an expansion known as the matricant, compare with Gantmacher [35]:

𝐅(𝑡,𝑡₀)=𝐈+^𝑡_𝑡0𝐋(𝜏)𝑑𝜏+^𝑡_𝑡0𝐋(𝜏)^𝜏_𝑡0𝐋(𝜎)𝑑𝜎𝑑𝜏+⋯. (6) The matricant has the following property:

𝐅(𝑡,𝑡₀)=𝐅(𝑡,𝑡₁)𝐅(𝑡₁,𝑡₀)(𝑡₀<𝑡₁<𝑡). (7) It follows that, if we split the interval [𝑡₀,𝑡] into 𝑛 parts by using intermediate points 𝑡₁,…,𝑡_𝑛−1 separated by Δ𝑡_𝑘=𝑡_𝑘−𝑡_𝑘−1 (𝑘=1,…,𝑛, with 𝑡_𝑛=𝑡),

𝐅(𝑡,𝑡₀)=𝐅(𝑡,𝑡_𝑛−1)⋯𝐅(𝑡₁,𝑡₀),with𝑡₀<𝑡₁<⋯<𝑡_𝑛−1<𝑡_𝑛=𝑡. (8) For an infinitesimally narrow time interval [𝑡_𝑘−1,𝑡_𝑘] we have by approximation

𝐅(𝑡_𝑘,𝑡_𝑘−1)=𝐈+𝐋(𝑡^∗_𝑘)Δ𝑡_𝑘+higherordertermsinΔ𝑡_𝑘, (9) in which 𝑡^∗_𝑘 is any point satisfying 𝑡_𝑘−1≤𝑡^∗_𝑘≤𝑡_𝑘. Equations (8) and (9) give rise to a representation in terms of a so-called multiplicative integral [35]:

𝐅(𝑡,𝑡₀)=∼^𝑡_𝑡0(𝐈+𝐋(𝜏)𝑑𝜏)_def=lim_{Δ𝑡𝑘→0}(𝐈+𝐋(𝑡^∗_𝑛)Δ𝑡_𝑛)⋯(𝐈+𝐋(𝑡^∗₁)Δ𝑡₁). (10) One recognizes the multiplicative counterpart of the Riemann sum approximation for ordinary (additive) integrals. One can show that this is identical to

𝐅(𝑡,𝑡₀)=^[35]∼^𝑡_𝑡0exp(𝐋(𝜏)𝑑𝜏)_def=lim_{Δ𝑡𝑘→0}exp(𝐋(𝑡^∗_𝑛)Δ𝑡_𝑛)⋯exp(𝐋(𝑡^∗₁)Δ𝑡₁). (11) Several properties of the deformation tensor are manifested in this representation. For instance, for square matrices 𝐀, 𝐁, one has

(i)det𝐀𝐁=det𝐀det𝐁, (ii)det(𝐈+𝜖𝐀)=1+𝜖tr𝐀+𝒪(𝜖²), (iii)detexp𝐀=exptr𝐀.

Consequently,

det𝐅(𝑡,𝑡₀)=∼^𝑡_𝑡0(1+tr𝐋(𝜏)𝑑𝜏)=∼^𝑡_𝑡0exp(tr𝐋(𝜏)𝑑𝜏). (12) In particular, a divergence free velocity field (div𝐯=tr𝐋=0) preserves volumes: det𝐅(𝑡,𝑡₀)=1. Furthermore, exp𝐀exp𝐁=exp(𝐀+𝐁) if [𝐀,𝐁]=0, whence for a stationary velocity field one obtains 𝐅(𝑡,𝑡₀)=exp((𝑡−𝑡₀)𝐋₀), as we already noticed. Finally, the multiplicative integral suggests a straightforward numerical approximation, namely by using either (10) or (11) without limiting procedure. In this case the two representations are of course no longer identical. For computational efficiency we have chosen the former, with Δ𝑡_𝑘=Δ𝑡 corresponding to the (constant) frame interval of our tagging MRI sequence, and 𝑡^∗_𝑘=𝑘Δ𝑡.

2.2. Strain

On the basis of the differential map 𝑓_∗, (2), and its transpose 𝑓^T_∗, one defines an intrinsic mapping on 𝑇𝑀_𝑋, known as the Lagrangian strain tensor [36],

1𝐸≡2(𝑓^T_∗∘𝑓_∗−id_𝑇𝑀𝑋)∶𝑇𝑀_𝑋→𝑇𝑀_𝑋∶𝑣↦𝐸(𝑣)=𝑣^𝑖𝐸^𝑗_𝑖𝐞_𝑗, (13) relative to a basis {𝐞_𝑗} of 𝑇𝑀_𝑋, with mixed tensor components

𝐸^𝑗_𝑖=12(𝑔^ℓ𝑗𝐹^𝛼_ℓℎ_𝛼𝛽𝐹^𝛽_𝑖−𝛿^𝑗_𝑖). (14) Here 𝑔^𝑖𝑗 are the components of the dual Euclidean metric tensor in domain coordinates {𝑋^𝑖,𝑀}, and ℎ_𝛼𝛽 are those pertaining to codomain coordinates {𝑥^𝛼,𝑁}. 𝐸 vanishes identically if 𝑓_∗ is an isometry, thus 𝐸 captures genuinely nonrigid deformations.

The general coordinate convention, see (14), which admits different bases for deformed and undeformed configurations, is instructive. Although one would normally prefer identical bases, let us consider the case in which the metric components ℎ_𝛼𝛽 are those induced from 𝑔_𝑖𝑗 by the deformation map itself (carry-along). That is, we assume that the coordinate frame for the deformed configuration is just the deformed coordinate frame of the reference configuration. One then expects the Lagrangian strain tensor to be nullified, because everything, including the local reference frames, is intrinsically deformed in a consistent manner. Indeed, if we write the infinitesimal material line element as

𝑑𝑠²=𝑔_𝑖𝑗𝑑𝑋^𝑖𝑑𝑋^𝑗=𝑔_𝑖𝑗𝜕𝑋^𝑖𝜕𝑥^𝛼𝜕𝑋^𝑗𝜕𝑥^𝛽𝑑𝑥^𝛼𝑑𝑥^𝛽≡ℎ_𝛼𝛽𝑑𝑥^𝛼𝑑𝑥^𝛽, (15) we recognize the deformation tensor, see (3), and it follows that

𝑔_𝑖𝑗=ℎ_𝛼𝛽𝐹^𝛼_𝑖𝐹^𝛽_𝑗. (16) As anticipated for the coordinate frames employed, (14) indeed reduces to

𝐸^𝑗_𝑖=12(𝑔^ℓ𝑗ℎ_𝛼𝛽𝐹^𝛼_ℓ𝐹^𝛽_𝑖(16)−𝛿^𝑗_𝑖1)=2(𝑔^ℓ𝑗𝑔_ℓ𝑖−𝛿^𝑗_𝑖)=0. (17) In practice one of course carries out computations relative to fixed (deformation independent) coordinate frames for undeformed and deformed configurations, but there is no reason to assume input and output frames to coincide a priori, nor to restrict oneself to Cartesian coordinate frames. Equation (14) permits us to use whatever bases we may consider convenient. This may be beneficial in certain situations, for example, those that suggest spherical, cylindrical or other coordinate systems depending on the configuration of acquisition data, compare with the phantom study by Young et al. [37], and the desired output representation.

Below we consider a single Cartesian coordinate system for both domain and codomain. (In such a system the representations of the various metric tensors, 𝑔_𝑖𝑗 and ℎ_𝛼𝛽 and their duals, 𝑔^𝑖𝑗 and ℎ^𝛼𝛽, all simplify to identity matrices.) Our analysis is confined to a single short-axis plane, whence 𝑛=2, that is, we only account for in-plane motion components. The theory trivially generalizes to 𝑛=3 and non-Cartesian coordinate grids.

3. Experiment

Figure 2 illustrates various scalar fields extracted from the strain tensor field, (13) and (14), by contraction with a pair of local unit vectors, namely radial (r = radial) and azimuthal (c = circumferential) basis vectors of the polar coordinate system centered at the midpoint of the region of interest, and those defining the strain tensor's eigensystem. If 𝑢,𝑣∈𝑇𝑀_𝑋 are two such unit vectors, then the local scalar quantity derived from the strain tensor is given by the inner product (𝑢,𝐸(𝑣))=(𝐸(𝑢),𝑣)=𝑔_𝑗𝑘𝑢^𝑘𝐸^𝑗_𝑖𝑣^𝑖. Tensor components are evaluated in Cartesian coordinates.

(a)

(b)

(c)

(d)

(e)

Figure 2 

Color-coded strain fields in short axis mid-slice cross-sections, time frames 1,6,11,16,21 (a full heart cycle subtends 42 frames), regularized through Gaussian convolution with spatial scale 𝜎=1.0, for one healthy volunteer. First row: Circumferential strain 𝐸cc. Second row: Radial strain 𝐸rr. Third row: Shear strain 𝐸cr. Fourth row: Minimal strain eigenvalue 𝐸min. Fifth row: Maximal strain eigenvalue 𝐸max.

Figure 3 illustrates temporal evolutions of these quantities spatially averaged over the respective regions of interest, together with their standard deviations. (Thus one should not confuse standard deviations with uncertainties in this plot.) Table 1 shows statistics for three healthy volunteers. Figure 4 shows that strain fields have the potential to detect local pathology.

Table 1 

Spatial average and standard deviation of strains from three healthy volunteers (labels I, II, III) over base, mid, and apex short-axis cross-sectional ROIs at that time frame at which 𝐸_𝑚𝑎𝑥 has attained its maximum in mid-slice (𝑡_𝐼=16, 𝑡_𝐼𝐼=15, 𝑡_𝐼𝐼𝐼=14). Rows correspond to, from top to bottom, 𝐸_cc, 𝐸_rr, 𝐸_cr, 𝐸_max, and 𝐸_min.

(a)

(b)

Figure 3 

Temporal evolution of scalar strain quantities over the first half of the heart cycle (i.e., mostly systolic) in mid-slice cross-section. Error bars indicate standard deviations over the spatial ROI in each time frame, thus capturing all sources of variation due to (primarily) actual spatial variability, noise, and numerics. Legends explain the various graphs. Notice the strong correlation between the extrinsic (polar system related) and instrinsic (eigensystem related) strains. (In the eigensystem, shear strain vanishes identically.)

(a)

(b)

Figure 4 

Radial (a) and circumferential (b) colour-coded strains in a basal short-axis slice at approximately end-systole for a patient. Red arrows on the left indicate regions where the respective strain values deviate from normal. These regions were confirmed to be static by a cardiologist. The patient's anamnesis mentions small infarcts on either side of these regions. Their locations, indicated by black arrows, are highlighted in late-enhancement MRI.

4. Conclusion and Future Work

We have proposed a novel, simple and robust linear model for extracting myocardial deformation and strain. Material motion and gradient velocity are determined by a multiscale linear system of algebraic equations for the phase images extracted from the tagging MRI data. These phase images are scalars, not densities, and are not hampered by tag fading due to 𝑇₁ relaxation.

We have analytically solved the linear matrix-ODE governing myocardial deformation. By discretizing the closed-form solution we have subsequently solved for the induced Lagrangian strain tensor field, yielding results that are typical for healthy volunteers, compare with Garot et al. [14], and atypical for a patient with a medical history of small infarcts on either side of a deviatory region and confirmed by late-enhancement MRI. This demonstrates the feasibility of our method. An advantage of our method is that only off-the-shelf algorithms are needed. This serves clarity, facilitates optimization, and enables implementation on dedicated hardware. The method is applicable in any spatial dimension, does not require conversion to a Cartesian sampling grid, is optimally conditioned due to local scale selection, and is readily adapted to other modalities such as velocity encoding MRI.

Our quantitative results in this feasibility study demand further evaluation so as to reveal their statistical significance and clinical value. Furthermore, experiments based on carefully controlled synthetic or phantom data will allow us to disclose the physical significance and numerical tolerance of the mathematical framework in terms of ground truth deformation properties. Moreover, by virtue of the transparent mathematical theory, it is also possible to assess tolerances on the basis of theoretical error propagation models [38], which can then be compared to those established experimentally for synthetic and phantom data. Consistency will greatly increase the confidence of the method for in vivo studies. Finally, although the local scale selection criterion exploited in our approach guarantees optimal conditioning, it does not necessarily yield optimal results due to the intrinsically parallel and thus spatiotemporally uncorrelated nature of the selected neighbouring scales. The effect of spatiotemporal regularization on the pattern of locally selected scales needs to be investigated. Future investigation will be needed to pursue all these recommendations.

Appendix

Data Acquisition Details

The short-axis MR tagging image data used in this study were acquired with a Philips Intera 1.5T scanner from three volunteers (Table 1) and a patient in basal, midventricular, and apical slices. The patient (Figure 4) had a history of severe stenoses, and small infarction areas confirmed with late-enhancement MRI 14 minutes after gadolinium contrast injection. A 2D multishot gradient-echo with Echo Planar Imaging (EPI factor 9) with breath-holding in end-expiration was used. Scan parameters were: TE 4.4ms, TR 19ms, flip angle 10^∘, field-of-view 300mm, scan matrix 128, acquisition voxel size 2.34×2.68×8mm³ reconstructed into 1.17×1.17×8mm³. Tagline spacing was 8mm.

Acknowledgments

The Netherlands Organisation for Scientific Research (NWO) is gratefully acknowledged for financial support. Jos Westenberg has provided the MRI data used in our experiments.