Multi-frequency Rayleigh damped elastography: in silico studies

G Model JJBE-2573; No. of Pages 13

ARTICLE IN PRESS Medical Engineering & Physics xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy

Multi-frequency Rayleigh damped elastography: in silico studies Andrii Y. Petrov ∗ , Paul D. Docherty, Mathieu Sellier, J. Geoffrey Chase Department of Mechanical Engineering, University of Canterbury, Christchurch, New Zealand

a r t i c l e

i n f o

Article history: Received 12 February 2014 Received in revised form 12 September 2014 Accepted 18 October 2014 Keywords: Magnetic resonance elastography Rayleigh damping Simulation studies Simultaneous multi-frequency inversion Model identifiability Mechanical properties

a b s t r a c t Rayleigh damping (RD) is commonly used to model energy attenuation for analyses of structures subjected to dynamic loads. In time-harmonic Magnetic Resonance Elastography (MRE), the RD model was shown to be non-identifiable at a single frequency data due to the ill-posed nature of the imaginary components describing energy dissipation arising from elastic and inertial forces. Thus, parametrisation or multi-frequency (MF) input data is required to overcome the fundamental identifiability issue of the model. While parametrisation allows improved accuracy of the identified parameters, simultaneous inversion using MF input data is a prerequisite for theoretical identifiably of the model. Furthermore, to establish good practical identifiability, frequencies should be separated over a wide range to produce different dynamic response. This research investigates the effects on practical identifiability of the RD model using MF data over different combinations of frequencies in noise-free heterogenous simulated geometry and compares the outcomes to reconstruction result based on single frequency input data. We tested eight frequencies in silico for a phantom type geometry comprises three independent material regions characterised by different mechanical properties. Combinations of two near or well separated frequencies are used to test the separation necessary to obtain accurate results, while the use of four or eight simultaneous frequencies is used to assess robustness. Results confirm expected non-identifiability of the RD model given single frequency input data. Practical identifiability of the RD parameters improved as more input frequencies were used for simultaneous inversion and when two frequencies were well separated. Best quality reconstruction results were achieved using full range data comprising eight available frequencies over a wide range. The main outcome is that high quality motion data over at least two frequencies over a wide range is required for establishing minimal practical identifiability of the model, while quality of the practical identifiability increases proportionally with more input frequencies used. Further simulation studies are required to determine acceptable signal-to-noise ratio (SNR) thresholds in motion data for accurate inversion of the RD parameters. © 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

1. Introduction Magnetic Resonance Elastography (MRE) quantitatively measures induced mechanical waves via phase-contrast MR imaging techniques [1,2]. Depending on the inverse problem methodology and constitutive model used, a number of parameters, describing mechanical properties of the material, can be estimated from the measured response data. Different models allow different parameters to be identified. For example, traditional linearly elastic [3–5]

∗ Corresponding author. Tel.: +1 8086700520. E-mail addresses: [email protected], [email protected] (A.Y. Petrov), [email protected] (P.D. Docherty), [email protected] (M. Sellier), [email protected] (J.G. Chase).

models can produce locally resolved maps of elasticity, also known as the stiffness or storage modulus (R ), while viscoelastic (VE) models [6–8] can provide additional information about viscosity or loss modulus (I ) of the material. However, more complex models, such as bi-phasic poroelastic [9] and Rayleigh damping (RD) models, incorporate additional parameters and, therefore, require higher quality motion data with adequate signal-to-noise ratio (SNR), as well as robust inversion methods for accurate parameter identification. Model fit to the observed response is critical to accurate parameter identification. However, true parameter estimation cannot be assured only by data-model match condition. Structural model identifiability is another absolute prerequisite to ensure unique and correct parameter estimation [10,11]. While VE models can produce accurate estimation of the parameters given a single

http://dx.doi.org/10.1016/j.medengphy.2014.10.007 1350-4533/© 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13 2

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

frequency of response data, RD models describe frequency dependent damping behaviour and rely on multi-frequency (MF) input data. Thus, accurate identification of the RD parameters is not possible given a single frequency of response data [12]. One way to overcome this limitation is to impose significant a priori information with a number of optimisation constraints and regularisation techniques under certain conditions and assumptions [13]. For example, parametrisation can be employed when only single frequency data is available [14]. However, this can only guarantee a limited trade off between data consistency and accuracy of estimated parameters. In fact, parametric approach is only recommended when a general qualitative profile of VE properties is sufficient. Damping has long been recognised as an important property in characterisation of the tissue structure and composition. It is directly associated with the viscosity of the material and has been recently linked to various pathologies in different organs, such as the liver [15,16], the breast [17,18], muscle tissues [19,20], and the brain [8,21–25]. Depending on micro structural configuration tissue would undergo different behaviour under applied mechanical stress. For example, fluidic or tightly arranged cells will display different attenuation behaviour due to changes in physiological consistency and vasculature, such as may arise with cancerous tumours [26]. Thus, ability to accurately capture and quantify damping information within the tissue of interest might provide useful diagnostic merit and help in early disease detection and investigation of the pathological progression over time. RD models are commonly used in structural dynamics, including automotive and civil structures [27]. They effectively extend the VE model, where damping is only defined by elastic forces, with the damping matrix represented as a scaled proportional combination of mass and stiffness terms. Thus, compared to the traditional VE model, the RD model incorporates an additional parameter to account for damping arising from viscous or inertial effects [28]. Previous phantom studies [29,30] showed sensitivity of the parameter to the material composition and ability to differentiate saturated porous tofu material from more tightly aligned course structure of the gelatine. Thus, it might provide valuable information for characterising tissue anatomy and might be particularly useful for imaging highly saturated structures, such as biological tissue. Since RD model renders two attenuation mechanisms governed by two different rheological effects, both elastic and inertial, it is hypothesised that the RD model can provide a more complete description of the intrinsic non-linear damping behaviour observed in the tissue in vivo. RD models have been previously investigated in silico [29] and in phantom studies [31]. The latter only evaluated identification of elasticity and overall damping, while accurate identification of imaginary components was not assessed. Hence, it was structurally identifiable at that level. Further phantom experiments [30] demonstrated challenges of the MF approach and confirmed the necessity of the wide separation of the input frequencies, as well as selection of the appropriate frequency dependent model to account for the dispersion characteristics of the given materials. In this context, preliminary simulation studies demonstrated the potential to correctly identify RD parameters on simulated noise-free motions at a single frequency within a homogeneous geometry with a high level of damping present [29]. However, more complex heterogeneous geometries with varying spatial damping profile governed by different attenuation mechanisms have not been properly investigated. This research determines the requirements for the identifiability of the RD model in heterogeneous simulated geometry using single frequency and MF noise-free data for image reconstruction processing.

2. Materials and methods 2.1. RD model in time-harmonic MRE The RD model is implemented through finite element (FE) based solution of a nearly incompressible linear isotropic Navier’s equation, defined as 2

∇ · ((∇ u + ∇ uT )) − ∇ (∇ · u) − ∇ P = −

∂ u , ∂t 2

(1)

where u is the complex column displacement vector within the medium;  is the first Lamé parameter ( = − 1/3 Pa for a nearly incompressible case),  is the second Lamé parameter, also known as a shear stiffness;  is the density of the material, ∇ P is a pressure term, related to volumetric changes through the bulk modulus, K, via the relationship: ∇ P = K ∇ · u. It can be discretised and defined as ¨ + Cu˙ + Ku = f Mu

(2)

for the mass, damping and stiffness matrices, M, C and K, respectively; displacement vector, u, and known sinusoidal input forcing, f. The RD assumption is facilitated through the RD definition: C = ˛M + ˇK.

(3)

For a time-harmonic case, where input and resulting response iωt , Eq. (2) in the frequency ˆ are f(x, t) = ˆf(x)eiωt and u(x, t) = u(x)e domain is written as





−ω2 1 −

i˛ ω



 



ˆ = ˆf. M + 1 + iωˇ K u

(4)

By assuming stiffness and density to be complex, so that ∗ = R + iI and ∗ = R + iI , Eq. (4) can be further simplified: ˆ = ˆf, [ − ω2 ∗ M + ∗ K ]u M

(5)

K

where = 1/(M) and = 1/(K) are normalised mass and stiffness matrices, respectively. In particular, R and R describe the real valued shear modulus and density in the undamped system. In contrast, I and I represent two different damping components related to elastic and inertial effects, respectively. They both can be expressed in terms of the RD parameters: I = ωˇR ,

I =

−˛R . ω

(6)

Rheologically, R and I represent storage and loss modulus of the material, while I is hypothesised to reflect fluid perfusion in a poroelastic media, such as biological tissue [12,30]. The resulting damping ratio,  d , is defined as d =

1 2



ˇω +

˛ ω



⇒

d =

1 2



I

R

−

I R



,

(7)

where the stiffness proportional term (ˇω) contributes damping linearly proportional to the response frequency and the mass proportional term (˛/ω) contributes damping inversely proportional to the response frequency. Therefore, the ˛ and ˇ coefficients govern the response to lower and higher input frequencies, ω, respectively. Previous studies found correlation between biological tissue response to various excitation frequencies in a form of power-law (PL). More specifically, Fabry et al. [32] observed PL dependency of the complex shear modulus in variety of individual cell types subjected to number of excitation frequencies over a wide range (10–300 Hz). Similarly, Robert et al. [33] studied ex vivo liver specimens using MF based MRE in the range between 40 and 100 Hz and found PL trend in the dynamic response of the complex shear modulus. A PL implies frequency dependant behaviour of the VE parameters, such as complex shear modulus, which can be derived from the fractional Zener and Kelvin VE models [34]. In this context, PL

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model

ARTICLE IN PRESS

JJBE-2573; No. of Pages 13

A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

3

derivation for the imaginary density can be assumed from the RD model formulation. Since the RD model incorporates an additional damping parameter that represents inertial attenuation and has a dominant effect in the lower frequency range, as well as inverse proportionality to the input frequency, the applied PL will consequently have inverse frequency relationship, as well. In this study, the PL relationship for all three reconstructed parameters (R , I and I ) is defined as R (ω) = R0 · ω ,

I (ω) = I0 · ω ,

I (ω) = I0 ·

 1  ω

,

(8)

where the static property values (at 0 Hz) R0 , I0 and I0 and a power component  are calculated conjointly for each node within the imaging volume. Collecting real and imaginary terms from Eq. (5), yields: ˆ = ˆf. [(R − ω2 R ) + i(I − ω2 I )]u

(9)

In Eq. (9) the force and displacement are complex valued, ˆf = ˆ = uR + iuI . Thus, by further collecting the real and fR + ifI and u imaginary terms, it can be rewritten as (R − ω2 R )uR − (I − ω2 I )uI = fR , 2

(10a)

2

(R − ω R )uI − (I − ω I )uR = fI .

(10b)

Multiplying Eqs. (10a) and (10b) by uR and uI , respectively, and summing them up, yields: (R − ω2 R )(uR · uR + uI · uI ) = fR · uR + fI · uI .

(11)

Given that the density (R ) can be assumed to be the same as density of water, Eq. (11) implies that R is uniquely identifiable as ˆ a direct function of the displacements: R = f [u]. In contrast, the shear wave attenuation of the material due to its damping effects is described by two parameters that have the same scaled model role. In particular, model behaviour defined by some particular values of I and I would be indistinguishable from I = I0 + x and I = I0 + x/ω2 . Hence, these parameters have a nonunique model role and cannot be uniquely identified with data from a single frequency without further a priori information [12]. 2.2. Elastographic reconstruction Reconstruction processing in this study computed by finite element (FE) based subspace non-linear inversion (SNLI) algorithm [35,36]. The algorithm operates on small overlapping subzones processed in a hierarchical order determined by the greatest error reduction in the objective function ( ), defined as (Â) =

Nm 

H

(um − uci (Â))(um − uci (Â)) , i i

(12)

Fig. 1. Simulated heterogenous geometry comprises 11 slices with two cylindrical inclusions governed by two different attenuation mechanisms, fully elastic and fully inertial. (a) For illustration purposes only three slices (top, middle and bottom), depicted in green coloured nodes, are shown, while all slices of both inclusions depicted in blue coloured nodes are superpositioned in the figure. (b) Full-field timeharmonic displacements at 1 mm amplitude, depicted in red coloured nodes, are applied at 25 Hz in the direction perpendicular to the orientation of the inclusions with other directions left unconstrained. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

i=1

where um represents the complex valued measured displacei

ment data at ith measurement point, uci (Â) is the co-located displacement calculated by a forward simulation of the model using a current estimate of the properties (Â), Nm is the number of measurements and the superscript H denotes the complex conjugate transpose. The conjugate-gradient (CG) method is used to update the material property distribution to improve agreement between FE computed displacements and the displacements captured by MRI. MF inversion performs simultaneous identification on motion data from two or more input frequencies. In the MF mode, the objective function ( MF ) is computed as a sum across all input frequencies and measurement points: MF (Â) =

Nω i Nm   ωi =1 j=1

H

{(um (ωi ) − ucj (ωi , Â))(um (ωi ) − ucj (ωi , Â)) }, (13) j j

where um (ωi ) represents the complex valued measured displacej ment data at ith frequency and jth measurement point, ωi , ucj (ωi , Â) is the analogous displacement computed by a forward computational model based on the current estimate of the properties (Â) at that particular frequency. The material properties vector  across two or more frequencies is computed using a frequency dependant PL model.

2.3. In silico studies Simulation studies were performed to investigate the ability of the SNLI reconstruction algorithm using the RD material model to correctly identify these three mechanical properties on a heterogeneous geometry. The geometry had a rectangular shape with two embedded cylindrical inclusions of the same size (Fig. 1(a)). The FE mesh for forward calculations consisted of 20 097 nodes comprising 2167 elements encompassed in 11 slices. The computed

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model

ARTICLE IN PRESS

JJBE-2573; No. of Pages 13

A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

4

Table 1 The material property distribution in the simulated heterogenous geometry. Parameter Units

R (Pa)

I (Pa)

I (kg/m3 )

d (%/100)

Background Inclusion 1 Inclusion 2

3300 3300 3300

165 165 1650

−50 −500 −50

5 27.5 27.5

displacement were supported on 27 noded quadratic hexahedral elements, while the material properties were supported on eight noded linear hexahedral elements with 1 mm isotropic resolution. To allow higher damping levels within the inclusions, compared to the background, one inclusion (inclusion 1 in the lower left corner) had I values increased by the factor of 10, and the other (inclusion 2 in the upper right corner) had I values increased by 10×. Thus, inclusion 1 was inertially damped, while inclusion 2 was elastically damped. These increased values in both imaginary components within the inclusions corresponded to damping levels of 27.5%, compared to 5% within the background. The material property distribution for the background and inclusions were specified as given in Table 1. Time-harmonic excitation at multiple frequencies (5, 10, 20, 25, 30, 40, 60 and 70 Hz) was applied to this geometry in the direction perpendicular to the orientation of the inclusions with other directions left unconstrained. The resulting displacements were measured at nodal resolution. An example of the direction and amplitudes of the applied displacements to the mesh at 25 Hz is depicted in Fig. 1(b). The element edge length was approximately 1.5 mm, corresponding to 16 nodes per shear wavelength for the prescribed properties. No noise was added to the simulated motion data. Full three-parameter based RD reconstructions were performed with different combinations of frequencies to investigate the identifiability of the RD based elastographic reconstruction. The frequency combinations utilised were: single frequency (40 or 60 Hz), two combined frequencies over a narrow range (5 and 10 Hz) and two combined frequencies over a wide range (10 and 60 Hz). Furthermore, identifiability of the model was also assessed against four and eight frequencies of input. To allow a fair comparison between qualitative and quantitative trends in parameters behaviour across different combinations of multiple input frequencies, the reconstruction results were performed at 40 Hz. No regularisation methods were applied during MF inversion processing due to the well-posed nature of the imaginary components when MF inversion is applied. No spatial filtering (SF) was applied to single frequency inversion, while initial and final SF weights for remaining MF inversions were 0.25 and 0.15. The expected theoretical outcome of the RD reconstruction using MF inversion and noise-free simulated data was a clear delineation of the RD parameters. Since no noise was added to the simulated motion data, MF based RD reconstruction was expected to produce accurate qualitative and quantitative estimates for all three parameters (R , I , I ). Utilising eight frequencies over a wide range for reconstruction processing was expected to produce the most accurate high quality reconstructions with minimal error of estimated values. Using four different frequencies (two in the lower range and two in the higher range) was also tested to observe whether increased data availability improves delineation of the imaginary model parameters. MF inversion based on two input frequencies over a wide range is a prerequisite for establishing practical identifiability and thus, was also predicted to give marginally accurate RD reconstruction. However, the use of two frequencies over a narrow range was expected to produce reconstruction artefacts and inaccurate qualitative and quantitative estimates. To allow comparison between reconstruction results using MF inversion, all power-law based reconstructions were performed for 40 Hz.

Two single frequency cases at 40 and 60 Hz were also considered to confirm predicted outcomes of the previously documented structural identifiability analysis [12]. Total non-identifiability of the imaginary components was expected for single frequency based RD reconstruction, while locally resolved elasticity maps were anticipated to be recovered reasonably well, as predicted by Eq. (11). 2.3.1. Reconstruction protocol The SNLI algorithm was used to estimate material properties using simulated displacement data. Reconstruction computations were carried out on the High Performance Computing (HPC) system Blue Fern P575 (University of Canterbury, Christchurch, New Zealand). A total of 32 processors were employed using a parallel Message Passing Interface (MPI). The average runtime for the reconstruction processing for a single frequency data was 5 h. Runtime increased with more frequencies used. For example, utilising all eight frequencies resulted in higher computational intensity and 72 h of total computing to produce 60 iterations. Simultaneous MF inversion of displacement data using the PL model was performed with the following parameters: an isotropic subzone size of 24 mm × 24 mm × 24 mm with the subzone overlap of 0.15% × 0.15% × 0.15%. Displacements were approximated on the mesh with 1.8 mm × 1.8 mm × 1.9 mm voxel resolution, providing approximately 16 nodes per wavelength for the FE forward problem. Each reconstructed parameter was interpolated at different resolution levels [37]. Initial parameter estimates were: R = 2000 Pa; I = 50 Pa and I = −5 kg/m3 . The real density (R ) and the bulk modulus (K) were constrained from inversion and were set to constants: R = 1000 kg/m3 and K = 13250 Pa. Convergence was achieved for all material properties across the reconstruction domain. Stabilisation of statistical indicators (median and interquartile range (IQR)) of parameter behaviour across all nodes was used to evaluate of the convergence for a particular parameter. The background and inclusions were then segmented and analysed to evaluate the presence of statistically significant differences among different regions of interest (ROIs) characterised by different mechanical properties. 3. Results Fig. 2 shows the simulated geometry where location of the two inclusions can be clearly seen. Fig. 3 shows full three-parameter RD reconstruction results using single frequency data at 40 and 60 Hz,

Fig. 2. Schematic 2D multi-slice representation of the simulated heterogeneous geometry, comprises 11 identical slices, where the top slices from left to right and downwards correspond to the top slices and downwards of the 3D geometry in Fig. 1. The geometry consists of the three independent ROIs defined by different combinations of material properties: background (dark grey colour), inclusion 1 (circle in lower left corner depicted in white colour) and inclusion 2 (circle in the upper right corner depicted in the light grey colour).

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

5

Fig. 3. Image results for the three-parameter RD reconstruction based on single frequency inversion at 40 Hz (left column) and 60 Hz (right column): (a and b) storage modulus R image (Pa); (b and c) loss modulus I image (Pa); (d and e) imaginary density I image (kg/m3 ); and (g and h) damping ratio  d image (%/100).

while Fig. 4 shows median and IQR of the identified parameters for each slice throughout the simulated geometry. Here, while R and  d were reconstructed reasonably well, clear non-identifiability of the RD parameters, I and I , is observed, as expected.

Fig. 5 shows the three-parameter PL based RD reconstruction at 40 Hz using MF inversion over two input frequencies of 5 and 10 Hz, and 10 and 60 Hz. Fig. 6 shows the identified parameters across three ROIs throughout the imaging volume at each slice. The

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13 6

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

Fig. 4. Quantitative analysis of the estimated parameters in single frequency based three-parameter RD reconstructions at 40 HZ (left column) and 60 Hz (right column) for three different ROIs corresponding to the background (black line), inclusion 1 (red line) and inclusion 2 (blue line): (a and b) storage modulus R image (Pa); (c and d) loss modulus I image (Pa); (e and f) imaginary density I image (kg/m3 ), and (g and h) damping ratio  d image (%/100). Uncertainty bars represent the [5 25 50 75 95] percentiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

7

Fig. 5. Image results for the three-parameter RD reconstruction at 40 Hz using PL model based on MF inversion over two input frequencies of 5 and 10 Hz (left column), and 10 and 60 Hz (right column): (a and b) storage modulus R image (Pa); (b and c) loss modulus I image (Pa); (d and e) imaginary density I image (kg/m3 ); and (g and h) damping ratio  d image (%/100).

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13 8

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

Fig. 6. Quantitative analysis of the estimated parameters in MF based three-parameter RD reconstructions at 40 Hz using PL model over two input frequencies of 5 and 10 Hz (left column) and 10 and 60 Hz (right column) for three different ROIs corresponding to the background (black line), inclusion 1 (red line) and inclusion 2 (blue line): (a and b) storage modulus R image (Pa); (c and d) loss modulus I image (Pa); (e and f) imaginary density I image (kg/m3 ), and (g and h) damping ratio  d image (%/100). Uncertainty bars represent the [5 25 50 75 95] percentiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

9

Fig. 7. Image results for the three-parameter RD reconstruction at 40 Hz using PL model based on MF inversion over four input frequencies of 10, 15, 60 and 70 Hz (left column) and eight frequencies of 5, 10, 15, 20, 25, 40, 60, 70 Hz (right column): (a and b) storage modulus R image (Pa); (b and c) loss modulus I image (Pa); (d and e) imaginary density I image (kg/m3 ); and (g and h) damping ratio  d image (%/100).

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13 10

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

Fig. 8. Quantitative analysis of the estimated parameters in MF based three-parameter RD reconstructions at 40 Hz using PL model over four input frequencies of 10, 15, 60 and 70 Hz (left column) and eight input frequencies of 5, 10, 15, 20, 25, 40, 60, 70 Hz (right column) for the three different ROIs corresponding to the background (black line), inclusion 1 (red line) and inclusion 2 (blue line): (a and b) storage modulus R image (Pa); (c and d) loss modulus I image (Pa); (e and f) imaginary density I image (kg/m3 ), and (g and h) damping ratio  d image (%/100). Uncertainty bars represent the [5 25 50 75 95] percentiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model

ARTICLE IN PRESS

JJBE-2573; No. of Pages 13

A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

results revealed that the combination of two frequencies over a wide range (10 and 60 Hz) produced better qualitative, as well as quantitative trends for RD parameters, compared with the results computed using two input frequencies over a narrow range (5 and 10 Hz), where poor practical identifiability of the RD parameters was observed, matching expectations. Finally, Fig. 7 shows the three-parameter PL based RD reconstruction at 40 Hz using MF inversion over four input frequencies (10, 15, 60 and 70 Hz) and alternatively eight input frequencies (5, 10, 15, 20, 25, 40, 60 and 70 Hz). Clear improvement in the accuracy of the reconstruction results is observed, with RD parameters having near perfect identifiably. Fig. 8 shows the range of identified parameters. Table 2 summarises the quantitative results for the three material regions for all simulation studies.

4. Discussion 4.1. Single frequency based three-parameter RD reconstruction. Single frequency based full three-parameters RD reconstruction demonstrated a clear failure of the algorithm to accurately identify RD parameters, I and I (Fig. 3(c)–(f)), while R was reconstructed reasonably well (Fig. 3(a) and (b)) in both cases. Qualitatively, R reconstruction shows uniform distribution of the property values within the background with small variance, while the presence of the inclusions is noticeable. Quantitative analysis confirmed uniform parameter behaviour within the background across the reconstruction domain with moderate variance within the inclusions in two cases (Fig. 4(a) and (b)). Quantitative estimate values for R within the background was 3291 Pa, which is in good agreement with the true values of 3300 Pa (Table 2). However, quantitative estimates for inclusion 1 and inclusion 2 were 3024 Pa and 3582 Pa, respectively, for 40 Hz reconstruction and 3127 Pa and

11

3406 Pa, for 60 Hz reconstruction indicating inaccurate identification of quantitative estimates. While marginal identifiability of the RD parameters was achieved within the background with I values of 200 and 172 Pa for 40 and 60 Hz single frequency based reconstructions, there was a clear failure of the algorithm to accurately identify RD parameters within the inclusions. Sporadic trends in I behaviour were noted in these cases (Fig. 4(c) and (d)), while similar trends were noted in I behaviour (Fig. 4(e) and (f)). Interestingly, quantified analysis of I revealed a symmetric linear taper with the centre in the middle slice, which could be due to reduced damping information near vibration sources (Fig. 4(e) and (f)). This trend was observed in both single frequency identification cases. Finally, the reconstructed  d image was qualitatively accurate, correctly confirming increased damping ratio within the inclusions compared to the background (Fig. 4(g) and (h)). Smooth distribution of the property values was observed in a homogeneous background material. Quantitatively,  d values for the background ( d = 4.2%) were reasonably accurate with values for both inclusions of 21% and 18% for 40 Hz reconstruction and 20% and 19% for 60 Hz reconstruction, compared with true value of 27.5% (Table 2), which again indicates inaccurate estimation of damping values. 4.2. Two frequency based three-parameter RD reconstruction Significant improvement in the quality of the reconstruction results, compared to a single frequency based reconstructions, was observed using simultaneous inversion of two input frequencies (Figs. 5 and 6). This improvement is especially pronounced using two frequencies over a wide range, as predicted by previously reported structural analysis of the RD model [12]. These outcomes match the initial hypothesis.

Table 2 Quantitative ROI analysis of the reconstructed parameters for the simulated heterogenous geometry across single and multiple frequencies. Parameter Units

R (Pa) Median [IQR]

I (Pa) Median [IQR]

I (kg/m3 ) Median [IQR]

d (%/100) Median [IQR]

True values Background Inclusion 1 Inclusion 2

3300 3300 3300

165 165 1650

−50 −500 −50

5 27.5 27.5

40 Hz Background Inclusion 1 Inclusion 2

3291 [3235–3336] 3024 [2874–3142] 3582 [3368–3701]

200 [154–219] 797 [768–839] 910 [814–927]

−22 [−26 to −16] −171 [−183 to −138] −123 [−153 to −95]

4.2 [3.8–4.5] 21 [19.2–22.7] 18 [15.9–20.4]

60 Hz Background Inclusion 1 Inclusion 2

3285 [3230–3331] 3127 [3016–3263] 3406 [3297–3496]

172 [141–209] 769 [649–820] 821 [691–910]

−22 [−28 to −12] −160 [−178 to −127] −117 [−128 to −91]

4 [3.4–4.4] 20 [17–22] 19 [16–21]

MF (10 and 60 Hz) PL at 40 Hz Background Inclusion 1 Inclusion 2

3079 [3059–3125] 2924 [2890–2974] 3097 [3058–3136]

163 [123–222] 550 [498–620] 1043 [808–1192]

−51 [−75 to −27] −187 [−242 to −116] −79 [−98 to −57]

5.1 [4.6–6.1] 19 [15–22] 20 [16–23]

MF (5 and 10 Hz) PL at 40 Hz Background Inclusion 1 Inclusion 2

4383 [4332–4448] 4444 [4386–4478] 4822 [4687–4934]

136 [109–173] 599 [484–650] 2072 [1587–2286]

−137 [−147 to −125] −1124 [−1289 to −856] −569 [−737 to −385]

8.4 [8–8.7] 63 [48–70] 50 [35–58]

MF (10, 15, 50 and 60 Hz) PL at 40 Hz 3101 [3070–3145] Background 2935 [2888–2988] Inclusion 1 3190 [3079–3245] Inclusion 2

162 [120–221] 545 [484–624] 1095 [829–1261]

−52 [−75 to −26] −215 [−269 to −128] −85 [−115 to −55]

5.2 [4.7–5.9] 20 [16–23] 22 [16–24]

F (8 freq) PL at 40 Hz Background Inclusion 1 Inclusion 2

135 [102–177] 323 [293–346] 1025 [730–1172]

−54 [−70 to −38] −274 [−331 to −205] −92 [−110 to −74]

5 [4.7–5.4] 19 [15–22] 20 [15–24]

3054 [3020–3097] 2996 [2964–3045] 3113 [3080–3143]

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13 12

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

Using two input frequencies over a narrow range (5 and 10 Hz) resulted in marginally accurate qualitative trends in all parameters. However, quantitative estimates for all parameters were not accurate. Qualitatively, correct trends in reconstructed I and I images was observed with the I image correctly capturing higher values within inclusion 2, while the I image correctly captured higher values in inclusion 1. However, both I and I images still displayed the presence of the second inclusion which was not supposed to be there (Fig. 5(c) and (e)). These trends confirm the hypothesised marginal practical identifiability using two input frequencies over a narrow range. Simultaneous inversion using two input frequencies over a wide range (10 and 60 Hz) resulted in better quality RD reconstructions compared with two input frequencies over a narrow range. Reconstructed R , I and  d images were accurate, correctly depicting overall qualitative trends in parameter behaviour (Fig. 6(b), (f) and (h)). However, the I image unexpectedly captured the presence of the second inclusion by increased property values (Fig. 6(f)). Overall, a relatively small level of variance was present in the reconstructed parameters. These results confirm the original hypothesis and provide evidence that at least two frequencies over a wide range are required for establishing practical identifiability. This outcome results from the greater differentiation in measured response, real and imaginary, when the frequencies are more widely separated. However, it also evident that two frequencies might still be not precise enough for clinical diagnosis of biological tissues. 4.3. Four and eight frequency based three-parameter RD reconstruction Full three-parameter RD reconstructions using four and eight input frequencies showed the best quality results (Figs. 7 and 8). Qualitatively, all four parameters were reconstructed correctly, with R and  d having highest quality image and I and I having slightly lower quality. Using eight input frequencies over a wide range significantly improves the quality of the reconstruction results. Reconstructed I images correctly delineate the elastically damped inclusion in the right upper corner (Fig. 7(d)) while the I images mainly capture the inertially damped inclusion located in the left corner of the sample geometry (Fig. 7(f)). Compared to MF inversion based on four input frequencies (two frequencies in the lower range and two frequencies in the upper range), the eight frequency case resulted in better practical identifiability of the model. Low levels of variance were generally observed in all four parameters within the three ROIs with minimal boundary effects present. Quantitatively, both the four and eight frequency based reconstructions did not result in exact estimates for I , I and  d values (Fig. 8(b), (d) and (f)). However, they were accurate enough and clearly located the inclusions to enable ready diagnosis in a clinical setting where differences in values matter more than exact precision. 4.4. Limitations While near-perfect practical identifiably was achieved using eight input frequencies, the practical application of the MF based MRE remains challenging. Multiple excitation frequencies would be required to obtain the data necessary for MF inversion, which is usually difficult in a hospital setting considering the time involved. However, the results clearly highlight the potential capability within this use of the algorithm. Furthermore, if only two or more input frequencies are used, unique identification of RD parameters is possible. However, the two frequencies have to be widely enough separated to assure the

phase and magnitude of responses are different enough to permit practical identifiability. In this context, if chosen frequencies give similar dynamic response in phase and magnitude, the unique identification of the RD parameters is not guaranteed. In this study, the loss of exact quantifiable parameter values with eight input frequencies is thus likely due to not having widely enough separation of input frequencies. The range used here was limited only because the 10–60 Hz range is what is clinically practical. Hence, identification quality will depend on the data available and the ability of the RD model to capture the observed dynamics. The SNLI algorithm is also a limitation. Theoretically, in noise free cases, exact reconstruction should be possible when structural identifiability is assured, as in the MF cases. The differences from this idea may be attributed to the subzone optimisation approach used. However, the required supercomputing for even these cases shows that a global optimisation approach is not likely to be feasible even if it removed this issue. Additionally, power-law function might not be appropriate to allow accurate recovery of the frequency dependent properties of the parameters. 5. Conclusion This research investigated effectiveness of the multiple frequencies on the identifiability of the RD model in the context of simulation studies, as well as assessed both qualitative and quantitative accuracy of the estimated parameters. Simulated fullfield time-harmonic noise-free motions at different frequencies within a heterogeneous geometry were used for reconstruction processing. Overall, sporadic trends in the RD parameters behaviour confirmed predicted non-identifiability of RD model using single frequency data. Furthermore, it was generally noted that overall identifiability quality of the three-parameter RD reconstruction increased proportionally with the number of frequencies used for the simultaneous MF inversion. Simultaneous inversion based on two input frequencies over a wide range yielded marginal identifiably of the model, while two frequencies over a narrow range resulted in a poor quality of the reconstruction results. Using four frequencies over a broad range resulted in a good practical identifiability of the model, while highest quality of the three-parameter RD reconstruction was achieved using eight input frequencies over a broad range. Role of the funding source No related funding. Ethical approval No human experimentation was undertaken. Conflict of interests The authors have no conflict of interest. All authors agree to submission. References [1] Lewa CJ. Magnetic resonance imaging in the presence of mechanical waves: NMR frequency modulation, mechanical waves as NMR factor, local temperature variations. Spectrosc Lett 1991;24(1):55–67. [2] Muthupillai R, Lomas DJ, Rossman PJ, Greenleaf JF, Manduca A, Ehman RL. Magnetic resonance elastography by direct visualization of propagating acoustic strain waves. Science 1995;269(5232):1854. [3] Walsh EK, Schettini A. Elastic behavior of brain tissue in vivo. Am J Physiol 1976;230(4):1058. [4] Walsh E, Schettini A. Calculation of brain elastic parameters in vivo. Am J Physiol Regul Integr Comp Physiol 1984;247(4):R693–700.

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007

G Model JJBE-2573; No. of Pages 13

ARTICLE IN PRESS A.Y. Petrov et al. / Medical Engineering & Physics xxx (2014) xxx–xxx

[5] Ferrant M, Warfield S, Nabavi A, Jolesz F, Kikinis R. Registration of 3D intraoperative MR images of the brain using a finite element biomechanical model. In: Medical image computing and computer-assisted intervention-MICCAI 2000. Springer; 2000. p. 249–58. [6] Bilston LE, Liu Z, Phan-Thien N. Linear viscoelastic properties of bovine brain tissue in shear. Biorheology 1997;34(6):377. [7] Wilcox RK, Bilston LE, Barton DC, Hall RM. Mathematical model for the viscoelastic properties of dura mater. J Orthop Sci 2003;8(3):432–4. [8] Sack I, Beierbach B, Hamhaber U, Klatt D, Braun J. Non-invasive measurement of brain viscoelasticity using magnetic resonance elastography. NMR Biomed 2007;21(3):265–71. [9] Perriez P, Kennedy F, Van Houten E, Weaver J, Paulsen K. Magnetic resonance poroelastography: an algorithm for estimating the mechanical properties of fluid-saturated soft tissues. IEEE Trans Med Imaging 2010;29(3):746–55. [10] Pia Saccomani M, Audoly S, D’Angiò L. Parameter identifiability of nonlinear systems: the role of initial conditions. Automatica 2003;39(4):619–32. [11] Audoly S, Bellu G, D’Angio L, Saccomani M, Cobelli C. Global identifiability of nonlinear models of biological systems. IEEE Trans Biomed Eng 2001;48(1):55–65. [12] Petrov AY, Sellier M, Docherty PD, Chase JG. Non-identifiablity of the Rayleigh damping model in magnetic resonance elastography. Math Biosci 2013;246(November (1)):191–201. [13] Van Houten EE. Parameter identification in a generalized time-harmonic Rayleigh damping model for elastography. PLOS ONE 2014;9(4):e93080. [14] Petrov AY, Sellier M, Docherty PD, Chase JG. Parametric-based brain magnetic resonance elastography using a Rayleigh damping material model. Comput Methods Programs Biomed 2014;116:328–39. [15] Klatt D, Asbach P, Rump J, Papazoglou S, Somasundaram R, Modrow J, et al. In vivo determination of hepatic stiffness using steady-state free precession magnetic resonance elastography. Invest Radiol 2006;41(12):841. [16] Asbach P, Klatt D, Hamhaber U, Braun J, Somasundaram R, Hamm B, et al. Assessment of liver viscoelasticity using multifrequency MR elastography. Magn Reson Med 2008;60(2):373–9. [17] Sinkus R, Tanter M, Catheline S, Lorenzen J, Kuhl C, Sondermann E, et al. Imaging anisotropic and viscous properties of breast tissue by magnetic resonanceelastography. Magn Reson Med 2005;53(2):372–87. [18] Yin M, Talwalkar JA, Glaser KJ, Manduca A, Grimm RC, Rossman PJ, et al. Assessment of hepatic fibrosis with magnetic resonance elastography. Clin Gastroenterol Hepatol 2007;5(10):1207–13. [19] Ringleb SI, Bensamoun SF, Chen Q, Manduca A, An K-N, Ehman RL. Applications of magnetic resonance elastography to healthy and pathologic skeletal muscle. J Magn Reson Imaging 2007;25(2):301–9. [20] Basford JR, Jenkyn TR, An K-N, Ehman RL, Heers G, Kaufman KR. Evaluation of healthy and diseased muscle with magnetic resonance elastography. Arch Phys Med Rehabil 2002;83(11):1530–6.

13

[21] Green MA, Bilston LE, Sinkus R. In vivo brain viscoelastic properties measured by magnetic resonance elastography. NMR Biomed 2008;21(7):755–64. [22] Sack I, Beierbach B, Wuerfel J, Klatt D, Hamhaber U, Papazoglou S, et al. The impact of aging and gender on brain viscoelasticity. Neuroimage 2009;46(3):652–7. [23] Wuerfel J, Paul F, Beierbach B, Hamhaber U, Klatt D, Papazoglou S, et al. MR-elastography reveals degradation of tissue integrity in multiple sclerosis. Neuroimage 2010;49(February (3)):2520–5. [24] Streitberger K-J, Wiener E, Hoffmann J, Freimann FB, Klatt D, Braun J, et al. In vivo viscoelastic properties of the brain in normal pressure hydrocephalus. NMR Biomed 2011;24(4):385–92. [25] Zhang J, Green M, Sinkus R, Bilston L. Viscoelastic properties of human cerebellum using magnetic resonance elastography. J Biomech 2011;44(10):1909–13. [26] Xu L, Lin Y, Han JC, Xi ZN, Shen H, Gao PY. Magnetic resonance elastography of brain tumors: preliminary results. Acta Radiol 2006;48(3):327–30. [27] Hall JF. Problems encountered from the use (or misuse) of Rayleigh damping. Earthq Eng Struct Dyn 2006;35(5):525–45. [28] Liu M, Gorman D. Formulation of Rayleigh damping and its extensions. Comput Struct 1995;57(2):277–85. [29] McGarry M, Van Houten E. Use of a Rayleigh damping model in elastography. Med Biol Eng Comput 2008;46(8):759–66. [30] Petrov AY, Docherty PD, Sellier M, Chase JG. Multi-frequency inversion in Rayleigh damped magnetic resonance elastography. Biomed Signal Process Control 2014;13:270–81. [31] Van Houten E, McGarry M, Perrinez P, Perreard I, Weaver J, Paulsen K. Subzone based magnetic resonance elastography using a Rayleigh damped material model. Med Phys 2011;38:1993. [32] Fabry B, Maksym GN, Butler JP, Glogauer M, Navajas D, Taback NA, et al. Time scale and other invariants of integrative mechanical behavior in living cells. Phys Rev E 2003;68(4):041914. [33] Robert B, Sinkus R, Bercoff J, Tanter M, Fink M. A novel fractal model to explain the rheology of liver tissue using MR-elastography. In: Proceedings of 14th annual meeting of ISMRM. 2006. p. 2560. [34] Holm S, Näsholm SP. Comparison of fractional wave equations for power law attenuation in ultrasound and elastography. Ultrasound Med Biol 2014;40(4):695–703. [35] Van Houten E, Paulsen K, Miga M, Kennedy F, Weaver J. An overlapping subzone technique for MR-based elastic property reconstruction. Magn Reson Med 1999;42(4):779–86. [36] Van Houten EEW, Miga MI, Weaver JB, Kennedy FE, Paulsen KD. Threedimensional subzone-based reconstruction algorithm for MR elastography. Magn Reson Med 2001;45(5):827–37. [37] McGarry M, Van Houten E, Johnson C, Georgiadis J, Sutton B, Weaver J, et al. Multiresolution MR elastography using nonlinear inversion. Med Phys 2012;39:6388.

Please cite this article in press as: Petrov AY, et al. Multi-frequency Rayleigh damped elastography: in silico studies. Med Eng Phys (2014), http://dx.doi.org/10.1016/j.medengphy.2014.10.007