Table of Contents Author Guidelines Submit a Manuscript
Computational and Mathematical Methods in Medicine
Volume 2016 (2016), Article ID 9343017, 12 pages
Research Article

Rescaled Local Interaction Simulation Approach for Shear Wave Propagation Modelling in Magnetic Resonance Elastography

Department of Robotics and Mechatronics, AGH University of Science and Technology, Al. Mickiewicza 30, 30-059 Krakow, Poland

Received 13 November 2015; Revised 16 December 2015; Accepted 17 December 2015

Academic Editor: Po-Hsiang Tsui

Copyright © 2016 Z. Hashemiyan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Properties of soft biological tissues are increasingly used in medical diagnosis to detect various abnormalities, for example, in liver fibrosis or breast tumors. It is well known that mechanical stiffness of human organs can be obtained from organ responses to shear stress waves through Magnetic Resonance Elastography. The Local Interaction Simulation Approach is proposed for effective modelling of shear wave propagation in soft tissues. The results are validated using experimental data from Magnetic Resonance Elastography. These results show the potential of the method for shear wave propagation modelling in soft tissues. The major advantage of the proposed approach is a significant reduction of computational effort.

1. Introduction

Mechanical properties of tissues are one of the most significant indicators used for detection of various abnormalities in medical diagnosis. Tumors and other pathologies often exhibit values of elastic moduli that are significantly different from healthy tissues. It is well known that none of the classical medical approaches, such as Computed Tomography (CT), Magnetic Resonance Imaging (MRI), and Ultrasonography (US), are able to detect mechanical properties of tissues that are diagnosed by palpation [1, 2]. Elastography is used extensively in diagnostic applications (e.g., liver fibrosis or breast tumors detection [39]) due to flexibility and noninvasiveness. Since abnormal tissues are often stiffer than the normal ones, medical diagnosis can be achieved. Although the method was developed in the late 1980s [1012] the major breakthrough came in the mid 1990s when a dynamic approach to elastography was proposed [13]. A Motion-Encoding Gradient (MEG) was introduced to a conventional MRI system leading to Magnetic Resonance Elastography (MRE) [1317].

Modelling in elastography relies on direct and inverse problems. The former relates to measurements of tissue responses to applied stresses. The latter is related to estimation of unknown mechanical properties from measured mechanical responses. Both problems are formulated using physical laws, which provide equations that relate biomechanical properties, such as shear modulus, Poisson’s ratio, viscosity, nonlinearity, and poroelasticity, to measured mechanical responses. Accurate models are required to predict displacement responses to different mechanical excitations to solve the inverse problem. For simple setups the equations that describe the direct problem have been solved analytically [18]. A similar approach used for irregular domains of elastically heterogeneous tissues is not possible in practice. Consequently, numerical simulations are used to ease this task. Modelling is used in MRE applications in order to create forward models that capture complex mechanisms of wave propagation in soft tissues. Previous studies in this field include various finite difference (FD) [1719] and Finite Element methods (FE) [15, 2024]. FE modelling has been used in previous studies for visualization of ultrasonic wave propagation [2531], elasticity reconstruction [21, 32], and shear wave propagation analysis in gelatin phantoms [3339].

The paper aims to develop a full three-dimensional (3D) model of shear wave propagation in a gelatin phantom for MRE applications. Some primary investigation has been performed for the bulk wave propagation model based on the Local Interaction Simulation Approach (LISA) [40]. In contrast to the previous work, current investigation focuses on the guided wave propagation with rescaling procedure. The major novelty of the presented work relates to the application of the Local Interaction Simulation Approach (LISA) for guided wave propagation and a rescaling procedure for the LISA is proposed for shear wave propagation modelling. This major novelty is considered to tackle numerical problems.

Then the LISA model is developed to examine density, shear modulus, and shear wavelength in a gelatin phantom. This study proposes the rescaling solution method in order to avoid numerical problems, especially related to wave amplitude. Numerical simulation results are compared with FE simulation results and MRE experimental measurements from a soft tissue mimicking an agarose gelatin phantom.

2. Theoretical Background

Elastic wave propagation in an isotropic linear medium is governed by the momentum balance given aswhere is the divergence of stress tensor, is an external volume force, and represents particle acceleration vector. The constitutive equation that relates stresses to strains in a linear elastic solid is given aswhere is the Kronecker delta, represents material dilatation given by , and represent the Lamé constants for the material. The strain () is defined through the strain tensor using the following relationshipwhere represents particle displacement components. Combining (1)–(3) the equation of equilibrium, that is, [41],governs wave propagation in an infinite elastic space and for practical problems must be amended by appropriate boundary and initial conditions describing the problem. Boundary conditions increase the complexity of the problem since they give rise to the so-called guided wave propagation problem, where global wave propagation patterns, that is, modes, travel at different, and possibly frequency-dependent, speeds, as explained in [41]. It is well known that the solution to (4) can be found only for simple canonical problems. Numerical simulations are used for more complex scenarios.

3. Numerical Models

This section describes numerical models used for shear wave propagation in soft tissues. Firstly FE model was developed as a reference. Then a LISA model is described. The major focus is on a rescaling procedure that is used to avoid numerical discrepancies.

3.1. Finite Element Model

The FE model used in the current investigations was developed using the Marc Mentat 2013 software package. Following the work presented in [36], a 3D cylindrical container, with a diameter of 200 mm and thickness of 20 mm, was modelled using gel phantom material properties. The bottom of the cylinder was fixed in the direction (see Figure 1). Altogether 36 000 elements of 2 × 2 mm radial and axial element size and 200 element along the circle were used. The phantom was modelled as a homogenous isotropic elastic solid with Poisson’s ratio . Harmonic sinusoidal motion of 150 Hz was applied to the center of the top cylinder surface as an excitation. Three different elastic moduli () were investigated, that is, 30, 60, and 120 kPa, to study the relationship between shear wavelengths and shear moduli. Similarly, numerical simulations were performed using three different material density () values, , , and , for each Young’s modulus. Material damping was assumed to be zero.

Figure 1: Elementary discretization scheme used for wave propagation modelling in the LISA 3D [42].

The shear wavelength () in the FE model was obtained by estimating the distances between wave peaks directly from response waveforms to make a direct comparisons with the results presented in [36].

3.2. Local Interaction Simulation Approach Model
3.2.1. Background of the LISA Model

The LISA, previously used for wave propagation in complex media [4249], has been applied for MRE shear wave propagation modelling. The algorithm of the LISA model is based on an FD approximation of (4) which discretises any structure under investigation into a grid of cells. Similar discretization is also used in the time domain when modelling is performed. All material properties are assumed to be constant within each cell but may differ between cells. The algorithm can be derived from the elastodynamic wave equation [42]where , is the vector of particle displacements, is the stiffness matrix, and are the differential operators matrices for stress and strain, respectively, and is the density. A comma before the subscript in (5) denotes differentiation. The matrix contains stiffness components that depend on Young’s moduli and Poisson ratios. The structure is discretised into parallelepiped cells for the 3D LISA wave propagation simulation, as illustrated in Figure 1. The junction of the eight cells characterizes the nodal point . The second time derivatives across the eight cells are needed to converge towards a common value at the point . In order to calculate a spatial derivative in the eight surrounding cells to , the central difference scheme is utilized. Then to obtain the solution, stress continuity across adjacent cells is constrained.

The following iteration equations are acquired for each displacement component for a general orthotropic case [42, 49]where is grid spacing in th direction, is the time step, denotes the th column of the matrix of signs, and displacement components are taken at time and point () if not stated otherwise. The index in the summation formula conveys the sign of the summed component. A detailed derivation of LISA equations can be found in [42]. The local interaction nature of boundary conditions based on matching conditions in the LISA model is the major advantage over the FD-based methods, when used for wave propagation. The so-called Sharp Interface Model (SIM) is used to average physical properties at interface grid points, which represent intersections of four elementary cells. When wave propagation problems in complex media with complex boundaries are studied, the SIM provides more accurate results, as demonstrated in [4247].

Shear wave propagation in the 3D cylinder, already described in the previous section, was modelled using the LISA approach. Numerical simulations involved the same material properties, boundary conditions, and excitation frequencies as in the FE model described in Section 3.1. These parameters were set following previous investigations reported in [36]. The 3D cylinder was meshed using 1 × 1 × 1 mm elements. Altogether 16 003 000 elements were used in the LISA model.

3.2.2. Rescaling Procedure

When a numerical technique is used, such as LISA, for wave propagation simulation, various numerical errors must be accounted for. It is well known that for certain material parameter values elastic waves are quickly damped out making results interpretation cumbersome. This is illustrated in Figure 2(a), where the wave field along the radial direction of the phantom is shown at a single time instant. Clearly, elastic waves are quickly damped out and determination of the wavelength becomes difficult. This problem is a consequence of numerical model properties and can be solved when certain model parameters are modified to avoid numerical discrepancies. Numerical model properties can be investigated in details through the iteration equations analysis. First, the severity of numerical damping can be analysed by considering roots of characteristic polynomial of a numerical scheme at hand, directly related to the Courant-Friedrichs-Lewy stability condition [50]. The latter is frequently invoked in the context of wave propagation modelling as the model parameters are required to meet certain restrictions for the analysis to be stable. This concept can be also used to quantify scheme’s accuracy, as will be shown next.

Figure 2: (a) The original shear waveform exhibiting attenuated amplitudes for the density at frequency 150 Hz. (b) The amplification factor based on different density.

Soft tissues are highly demanding from computational point of view. From physical perspective it is well known that mainly transversally polarized waves propagate in these structures [41]. When numerical modelling is used both types of waves normally coexist. However, analyzing material properties characteristic to soft tissues (these properties are extraordinary when compared to solid media) the difference in longitudinal and shear wave velocities can be immediately found, reaching the ratio of 10. As a consequence, the shear wave component, which is of particular interest for MRE, propagates under conditions far from the stability limit. Namely, the roots associated with the characteristic polynomial drive the waves to decay.

This drawback can be resolved twofold: by reformulating constitutive relationships in order to eliminate the longitudinal wave component, or by manipulating model parameters to push the shear wave closer to the stability limit. In the following work the second approach was employed as this requires no intervention in the solver structure, maintaining the flexibility of the method to model wider class of materials (i.e., solid media and soft tissues).

The longitudinal and shear wave speeds can be expressedThese definitions show that the density is a parameter that uniformly influences both longitudinal and shear wave velocities. Hence, the approach presented in the paper aims at improving the model properties by rescaling wave speeds. Following the work presented in [36], scaled density is used in numerical simulations. This procedure can be explained using a 1D example of wave propagation. The major focus is on the stability and amplitude accuracy of LISA. The objective of this study is to obtain information about the effect of density on the LISA model. The 1D finite difference equation, involved in numerical simulations, can be expressed aswhere is the displacement, is the time step index, relates to the node position, is known as the (dimensionless) Courant (or Courant-Friedrichs-Lewy (CFL)) number, is the wave velocity, and , are related to the time step and the element size, respectively. The stability analysis by means of the Fourier transform is known as the von Neumann analysis [50]. This analysis allows for expression of a governing equation as a recurrence relation that is particularly useful for establishing stability conditions. The key idea is the analysis of the amplification polynomial of the scheme, which is obtained by applying the Fourier transform to the governing FD equation. Once the amplification polynomial is established, certain restrictions are put on its roots. Although the analysis presented is for a 1D case, the entire procedure can be easily extended to provide general stability conditions for higher dimensions.

When (10) is used and the stability condition is obtained for various parameters, stable and unstable conditions can be analysed for various values of density. It is important to note that it is beyond this paper to put all the equation and formula involved in this analysis. Potential readers are referred to [50] for further details. After obtaining the roots of amplification polynomial which are conjugate pairs of the same complex number, the magnitude is the same for both. The magnitude of the roots of amplification factor can be expressed aswhere is the amplification factor, is courant number, , and is wave number.

The amplification factor that governs the numerical damping is analysed in Figure 2(b). According to (10), the amplification factor formula, an analysis becomes unstable for , since each consecutive displacement value is increased. Perfect preservation of wave amplitude, that is, without numerical damping, occurs for . For real applications, for example, involving longitudinal and shear waves, is always less than 1 for at least one wave mode; hence waveforms are always numerically distorted. The results shown in Figure 2(b) indicate that when the density decreases, the amplification factor increases (consistently with (7) and (8)) and the instability can be reached in 1D finite difference model.

Numerical simulations of the gel phantom using the LISA approach were conducted and analysis was performed to investigate the effect of density scaling parameter on numerical stability for the 3D case. Various scaling values were selected and respective wave velocities calculated. The initial density was assumed as . Then, five different scaling parameters were selected as , , , , and . The influence of the scaling on model properties is summarised in Table 1. Therein, the limit velocity for the model can be calculated asCalculation of the amplification factor for a general 3D case is cumbersome; hence it is not provided in the table. However, the general conclusions can be inferred from the CFL numbers as described next.

Table 1: The effect of scaled density on wave propagation velocities and numerical stability.

The results in Table 1 show that by increasing the density wave velocities are reduced (see (7) and (8)). As a result the values of the Courant number are also reduced. Figure 3(a) shows LISA-based simulation results, that is, displacement patterns for different rescaled densities. (Where in the figure horizontal line corresponds to distance of wave propagation form center.)

Figure 3: The shear waveform patterns: (a) density-scaled models and (b) after inverse rescaling.

If the only influence on wave amplitude was due to amplification factor, amplitude drop should have been observed in the simulation. Previous work on the effect of the courant number on pulse distortion in 1D finite difference schemes [51] confirmed these observations. Interestingly, the displacement amplitude increases with the scaling parameter. The latter increase is related to the amount and rate of energy transfer through the excitation. Increased densities result in larger kinetic energies delivered even for unchanged excitations that are prescribed by displacements. This explains the amplitude increase observed in Figure 3(a).

It is clear that once wave propagation is simulated with scaled densities, an inverse spatial scaling procedure should be applied to the results to retrieve proper responses. This is accomplished by an inverse scaling procedure employed for the space sensor waveforms. Again (7) and (8) were employed and space (wavelength) signals were multiplied by the square root of the relevant scaling factors. The results, shown in Figure 3(b), illustrate that the wavelength of the original signal is recovered after the rescaling procedure.

To illustrate the approach, dispersion curves for respective rescaled models were calculated and used to recover the original waveforms. In the following analysis, the mode is considered, as it is the dominant mode in this frequency range. In Figure 4, dispersion curves for three different scaling parameters are given ( in figures correspond by scaling factor). By applying the scaling factor to the original density, it affects the dispersion curve, respectively, which causes a certain change to the wave number of every dispersion curve. The rescaling parameter, square root of – based scaling factor, is also then obtained by analyzing dispersion curves plot between simulated original and scaled density by comparing the ratio of wave number of scaling density by wave number of original density and it proved the efficiency of propose method. To show an example, the waveform (density ) and one rescaled waveform () together with the corresponding dispersion curves are shown in Figure 5. The wavelength for the original waveform is equal to 21 mm whereas the wavelength for the rescaled waveform is equal to 12 mm. The original waveform (shown in Figure 5(a)) can be recovered from the rescaled waveform (Figure 5(c)) when the latter is multiplied (scaled back) by the square root of the scaling factor (i.e., square root of 3 in this case) and vice versa. The results are shown in Figure 6. The analysis of dispersion curves reinforces the condition that the guided not bulk wave theory should be used in the case investigated, as discussed further in Section 5.

Figure 4: The dispersion plot for the original density and three , , and scaled densities.
Figure 5: Numerical simulations of shear wave propagation: (a) original waveform (density ); (b) dispersion curves corresponding to the original waveform; (c) rescaled waveform (density ); (d) dispersion curves corresponding to the rescaled waveform.
Figure 6: Simulated shear waveforms: (a) after rescaling from to ; (b) after rescaling from to .

In summary, two interesting observations can be made after the analysis performed in this section. Firstly, the wave amplitude increases when density is rescaled towards larger values. Secondly, the inverse rescaling of waveforms allows one to reproduce accurately the original wavelengths.

4. Magnetic Resonance Elastography: Experimental Data

The MRE data from the experiments reported in [36] were used as a reference in the current investigations. The phantom used in the experiment was a 3D cylinder filled with 2% agarose gel. The geometry of the cylinder was as follows: diameter 150 mm and height 20 mm. The MRE tests were conducted using the 1.5 T General Electric Signa CT scanner. The phantom was placed in a head coil and an electromechanical driver was placed on the top surface of the phantom in order to generate shear waves corresponding to the excitation frequency of 150 Hz. The experimental setup used is shown in Figure 7.

Figure 7: The experimental setup used for MRE tests: I. phantom; II. applicator of the electromechanical driver; III. electromechanical driver; IV. head coil [36].

The propagation of elastic waves in the phantom was imaged with an MRE pulse sequence sensitive to motion in the horizontal direction. The shear wavelength was estimated manually by calculating the distances between the adjacent wave peaks. Also, the mean of shear wavelength was measured by averaging the wavelength over the four phase offsets. Subsequently, for isotropic elastic infinite solid, an estimate of the local shear modulus can be obtained from the local estimate of wavelength as [13]The shear wavelength for phantom estimated by MRE was 38.00 ± 2.12 mm at 150 Hz. This corresponds to the mean value of 28.5 kPa for the shear modulus . The shear modulus was also estimated using a dynamic multifrequency shear test with the DMA 2980 machine for polymer testing to obtain the value of 30 kPa. The density was estimated experimentally as .

5. Numerical Simulation Results

Numerical simulations of shear wave propagation in the phantom described in Section 4 were performed using FE and LISA models. The results are presented in this section.

5.1. Shear Wave Propagation in Soft Tissue

Simulated FE, LISA, and experimental MRE shear wave propagation patterns are presented in Figures 8(a), 8(b), and 8(c), respectively. The simulated results were obtained for Young’s modulus and the density . The density scaling was applied in LISA models to avoid numerical problems related to excessive wave attenuation. Subsequently, the rescaling procedure was used in postprocessing to recover proper waveforms. The results in Figure 8 show that the simulated and experimental wave patterns reveal the same wavelengths. Small differences between FE and LISA models can be attributed to different formulations of the FE and LISA equations used and differences in meshes.

Figure 8: Shear wave propagation patterns for the agarose gel phantom: (a) FE mode [40]; (b) LISA model; (c) MRE experiment [36].

Subsequently, the out-of-plane displacement component responses were acquired from the simulated (FE and LISA models before and after scaling) and experimental (MRE measurements) data. The results, presented in Figure 9, show good agreement between simulated and experimental displacements. It is also important to note that after scaling the amplitude of the LISA model is improved. Next, the shear wavelength was computed from the distance between two successive peaks (or valleys). The wavelengths were estimated as = 37.5 mm and = 37 mm for the simulated FE and LISA models, respectively. These results correspond quite well with the MRE-based experimental value of the wavelength = 38 mm. However, the computational effort of 10 seconds the LISA model compares favorably if compared with the 2640 seconds for the FE model.

Figure 9: Shear wave propagation, comparison of displacement waveforms for the FE model and LISA model before and after scaling and MRE measurements.

Following these investigations, simulated shear wavelengths, calculated for different values of elastic moduli and density, were compared with the relevant analytical values calculated from (7) and (8) for the bulk wave propagation problem. Four different elastic moduli, that is, 30, 60, 90, and 120 kPa, and three different densities, that is, , , and , were investigated. Figure 10 presents the results for the 150 Hz excitation frequency. Here, the three continuous solid, dashed, and dotted curves give the values of shear wavelengths calculated from (12) for infinite medium propagation model.

Figure 10: Comparison of shear wavelengths estimated from the FE and LISA models of bulk wave propagation with the relevant analytically estimated wavelengths for different elastic moduli and material densities (0.5, 1 and ) at frequency of 150 Hz [40].

Although the results are quite consistent for lower values of elastic moduli, significant discrepancies between numerically (FE and LISA) and analytically (bulk wave propagation solution) estimated results can be observed for higher values of elastic moduli (corresponding to larger wavelengths), particularly for lower densities. These discrepancies are further discussed in the next section.

5.2. Guided Wave Propagation in Soft Tissues

Equations (12) provide the relationships between excitation frequency, wavelength, and elastic constants for an infinite elastic space. Thus any estimation of wavelengths, as discussed in the previous sections, and consequently estimation of elastic properties based on these wavelengths is accurate only for the infinite space assumption. This assumption can be approximately fulfilled for the following two conditions: () wavelength estimates are made sufficiently far from the object’s boundaries; () wavelengths are small when compared with distances from the boundaries. Both conditions can be achieved when excitation frequencies are selected to obtain sufficiently short wavelengths. However, near the boundaries wavelength estimations will not be accurate, unless the effect of the interfaces is taken into account.

Figure 11 shows the through-thickness cross-section of the wave propagation displacement field for the analysed phantom model. The results, obtained for the 150 Hz excitation frequency, show that the displacement varies across the thickness of the phantom, from a finite value (top) to zero (bottom). This nonuniform displacement distribution indicates that the wave field is strongly affected by the (top and bottom) boundaries; as a result the (global) wavelength is different from the wavelength for the assumed theoretical infinite space case.

Figure 11: Example of shear wave propagation displacement profile in the phantom.

These numerically estimated values were compared in Figure 10 with the theoretical values obtained for the assumed infinite space. The results show that the numerical and analytical results start to diverge for wavelengths larger than . This corresponds to the thickness of the phantom. Thus the analysed wave propagation field in the phantom corresponds to the guided wave field rather than to the bulk wave field (assumed in (12)). Guided wave propagation involves partial waves, that is, waves propagating in an infinite space that interact with the (top and bottom) boundaries. These partial waves undergo multiple reflections and mode conversions forming global displacement patterns, that is, wave modes. Therefore, wavelength estimation in the analysed model should involve the relevant dispersion equations for guided waves rather than (7) and (8) for bulk waves. This problem can be solved semianalytically or numerically, as illustrated in [48]. A semianalytical approach, based on the LISA iteration equations, was used in the current paper for wavelength estimation. The vertical () displacement component at the bottom surface of the phantom was constrained. The results for various Young’s moduli and densities are presented in Figure 12. This time, the wavelengths estimated from the guided wave propagation model are compared with the relevant wavelengths estimated from the FE infinite space model (i.e., from (12)). When the results are analysed two distinct wavelength ranges can be distinguished in Figure 12. The semianalytical solution for guided wave propagation compares very well with the bulk wave model for wavelengths shorter than . In contrast, the results for the guided and unbounded media differ significantly for longer wavelengths. The results in Figures 11 and 12 indicate that the guided wave propagation model rather than the bulk wave propagation model (that was originally employed in [36]) should be used for wavelength estimation in the case investigated.

Figure 12: Comparison of shear wavelengths estimated from the semianalytical LISA guided wave model with the relevant analytically estimated wavelengths from bulk wave propagation model for different elastic moduli and material densities () at frequency of 150 Hz.

Since guided wave propagation is inevitably associated with wave interactions with boundaries, the effect of boundary conditions was investigated. Two different boundary conditions, namely, fixed and free ends, were examined. Altogether five different model scenarios, for both the FE and semianalytical LISA models, were analysed: () 20 mm thick phantom with bottom surface fixed in the -direction; () 20 mm thick phantom with the -component fixed; () 20 mm thick phantom with the free boundaries; () 40 mm thick phantom with the -direction fixed; and () 40 mm thick phantom with the free boundaries. The relevant calculations were performed for Young’s modulus and density . All numerically simulated results are presented in Figure 13 and compared with the semianalytical guided wave results. The comparison of the results obtained for the 20 mm thick phantom between the first three model scenarios investigated shows that the estimated wavelength increases when the boundary is fixed. Significantly different responses are obtained particularly for the -constrained direction. This can be attributed to the dominant shear wave propagating in the phantom. In other words analysed displacements are mainly in the -direction and thus the response is more sensitive to this type of boundary condition; thus, a substantial increase in the wavelength can be observed for the model with the -displacement component constrained, as expected.

Figure 13: Five different model boundary condition, for both the FE and semianalytical LISA models; () 20 mm fixed, () 20 mm fixed, () 20 mm free, () 40 mm fixed, () 40 mm free (for Young’s modulus and density ).

When the 40 mm thick phantom is analysed, the wavelength is larger for free conditions, if compared with the relevant 20 mm thick phantom in the semianalytical model. The value of wavelength is then further increased by the constraint in the -direction (comparison of scenarios () and (), in the 40 mm case) but an opposite trend was observed for the FE model when boundaries are changed from the free to the constrained in -direction. However, it is important to note the wavelength was calculated in this case from the peak distances and was less accurate than the semianalytical solution to the dispersion relation.

6. Conclusions

A 3D rescaled LISA model has been proposed for shear wave propagation analysis. Numerical simulations have been performed to analyze the shear wavelength, that is, the primary parameter characterizing shear modulus, in order to examine several factors that influence shear modulus estimation in homogenous phantoms.

The results show that rescaled LISA can be used very efficiently for shear wave propagation modelling in MRE investigations. Good results agreements have been achieved between the LISA-based, FE model, and experimental MRE measurements. The major advantage of the proposed rescaled LISA method is computational efficiency. Significant reductions of computational effort have been achieved when compared with the classical FE modelling approach. The computational time was reduced more than 260 times for the case investigated.

The results also demonstrate that shear wavelength estimated from the presented LISA and FE models are reasonably close to the theoretical calculations, for homogenous elastic cylindrical phantoms investigated, for shorter wavelengths (i.e, for lower Young’s moduli and high densities). In contrast, the solutions based on guided wave propagation are more accurate for longer wavelengths. Also by the rescaling procedure which is presented in this paper, the wave amplitude problems related to numerical errors in soft tissues modelling can be avoided. This analysis can serve as an indicator of interfacial conditions for complex wave propagation in biological tissues.

Conflict of Interests

There is no conflict of interests involved.


The work presented in this paper was supported by the Foundation for Polish Science under the research WELCOME Project no. 2010-3/2 (Innovative Economy, National Cohesion Programme funded by EU). The fourth author would like to acknowledge research financial support from the Faculty of Mechanical Engineering and Robotics, AGH University of Science and Technology. The authors would also like to thank Professor Kai-Nan An from Mayo Clinic College of Medicine in Rochester, USA, and Dr. Frank Chen from Excelen in Minneapolis, USA, for experimental MRE data used in the paper.


  1. Y. K. Mariappan, K. J. Glaser, and R. L. Ehman, “Magnetic resonance elastography: a review,” Clinical Anatomy, vol. 23, no. 5, pp. 497–511, 2010. View at Publisher · View at Google Scholar · View at Scopus
  2. S. Kathiravan and J. Kanakaraj, “A review on potential issues and challenges in MR imaging,” The Scientific World Journal, vol. 2013, Article ID 783715, 10 pages, 2013. View at Publisher · View at Google Scholar · View at Scopus
  3. M. Yin, J. A. Talwalkar, K. J. Glaser et al., “Assessment of hepatic fibrosis with magnetic resonance elastography,” Clinical Gastroenterology and Hepatology, vol. 5, no. 10, pp. 1207.e2–1213.e2, 2007. View at Publisher · View at Google Scholar · View at Scopus
  4. L. Huwart, C. Sempoux, N. Salameh et al., “Liver fibrosis: noninvasive assessment with MR elastography versus aspartate aminotransferase-to-platelet ratio index,” Radiology, vol. 245, no. 2, pp. 458–466, 2007. View at Publisher · View at Google Scholar · View at Scopus
  5. J. Chen, M. Yin, K. J. Glaser, J. A. Talwalkar, and R. L. Ehman, “MR elastography of liver disease: state of the art,” Applied Radiology, vol. 42, no. 4, pp. 5–12, 2013. View at Google Scholar · View at Scopus
  6. D. B. Plewes, J. Bishop, A. Samani, and J. Sciarretta, “Visualization and quantification of breast cancer biomechanical properties with magnetic resonance elastography,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1591–1610, 2000. View at Publisher · View at Google Scholar · View at Scopus
  7. A. L. Plewes, J. Mcknight, J. L. Kugel et al., “MR elastography of breast cancer: preliminary results,” American Journal of Roentgenology: Diagnostic Imaging and Related Sciences, vol. 178, no. 6, pp. 1411–1417, 2002. View at Google Scholar
  8. R. Sinkus, M. Tanter, S. Catheline et al., “Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography,” Magnetic Resonance in Medicine, vol. 53, no. 2, pp. 372–387, 2005. View at Publisher · View at Google Scholar · View at Scopus
  9. L. Q. T. Afroze, Developing Magnetic Resonance Elastography (Mre) Breast Actuation System for Detecting Breast Cancer, University of Canterbury. Mechanical Engineering, 2012.
  10. R. M. Lerner and K. J. Parker, “Soneoelasticity images in ultrasonic tissue characterization and echographic imaging,” in Proceedings of the 7th European Communities Workshop, J. Thijssen, Ed., Nijmegen, The Netherlands, 1987.
  11. R. M. Lerner, k. J. Parker, J. Holen, R. Gramiak, and R. C. Waag, “Sono-elasticity: medical elasticity images derived from ultrasound signals in mechanically vibrated targets,” in Acoustical Imaging, vol. 16 of Acoustical Imaging, pp. 317–327, Springer, New York, NY, USA, 1988. View at Publisher · View at Google Scholar
  12. J. Ophir, I. Céspedes, H. Ponnekanti, Y. Yazdi, and X. Li, “Elastography: a quantitative method for imaging the elasticity of biological tissues,” Ultrasonic Imaging, vol. 13, no. 2, pp. 111–134, 1991. View at Publisher · View at Google Scholar · View at Scopus
  13. R. Muthupillai, D. J. Lomas, P. J. Rossman, J. F. Greenleaf, A. Manduca, and R. L. Ehman, “Magnetic resonance elastography by direct visualization of propagating acoustic strain waves,” Science, vol. 269, no. 5232, pp. 1854–1857, 1995. View at Publisher · View at Google Scholar · View at Scopus
  14. A. Manduca, T. E. Oliphant, M. A. Dresner et al., “Magnetic resonance elastography: non-invasive mapping of tissue elasticity,” Medical Image Analysis, vol. 5, no. 4, pp. 237–254, 2001. View at Publisher · View at Google Scholar · View at Scopus
  15. J. Bishop, A. Samani, J. Sciarretta, and D. B. Plewes, “Two-dimensional MR elastography with linear inversion reconstruction: methodology and noise analysis,” Physics in Medicine and Biology, vol. 45, no. 8, pp. 2081–2091, 2000. View at Publisher · View at Google Scholar · View at Scopus
  16. J. B. Weaver, E. E. W. Van Houten, M. I. Miga, F. E. Kennedy, and K. D. Paulsen, “Magnetic resonance elastography using 3D gradient echo measurements of steady-state motion,” Medical Physics, vol. 28, no. 8, pp. 1620–1628, 2001. View at Publisher · View at Google Scholar · View at Scopus
  17. R. Sinkus, J. Lorenzen, D. Schrader, M. Lorenzen, M. Dargatz, and D. Holz, “High-resolution tensor MR elastography for breast tumour detection,” Physics in Medicine and Biology, vol. 45, no. 6, pp. 1649–1664, 2000. View at Publisher · View at Google Scholar · View at Scopus
  18. M. M. Doyley, “Model-based elastography: a survey of approaches to the inverse elasticity problem,” Physics in Medicine and Biology, vol. 57, no. 3, pp. R35–R73, 2012. View at Publisher · View at Google Scholar · View at Scopus
  19. K. R. Raghavan and A. E. Yagle, “Forward and inverse problems in elasticity imaging of soft tissues,” IEEE Transactions on Nuclear Science, vol. 41, no. 4, pp. 1639–1645, 1994. View at Publisher · View at Google Scholar · View at Scopus
  20. K. J. Parker, S. R. Huang, R. A. Musulin, and R. M. Lerner, “Tissue-response tomechanical vibrations for sonoelasticity imaging,” Ultrasound in Medicine & Biology, vol. 16, no. 3, pp. 241–246, 1990. View at Google Scholar
  21. E. E. W. Van Houten, M. I. Miga, J. B. Weaver, F. E. Kennedy, and K. D. Paulsen, “Three-dimensional subzone-based reconstruction algorithm for MR elastography,” Magnetic Resonance in Medicine, vol. 45, no. 5, pp. 827–837, 2001. View at Publisher · View at Google Scholar · View at Scopus
  22. A. Samani, J. Bishop, and D. B. Plewes, “A constrained modulus reconstruction technique for breast cancer assessment,” IEEE Transactions on Medical Imaging, vol. 20, no. 9, pp. 877–885, 2001. View at Publisher · View at Google Scholar · View at Scopus
  23. M. I. Miga, “A new approach to elastography using mutual information and finite elements,” Physics in Medicine and Biology, vol. 48, no. 4, pp. 467–480, 2003. View at Publisher · View at Google Scholar · View at Scopus
  24. J. C. Brigham, W. Aquino, F. G. Mitri, J. F. Greenleaf, and M. Fatemi, “Inverse estimation of viscoelastic material properties for solids immersed in fluids using vibroacoustic techniques,” Journal of Applied Physics, vol. 101, no. 2, Article ID 023509, 2007. View at Publisher · View at Google Scholar · View at Scopus
  25. S. R. Ghorayeb, E. Maione, and V. La Magna, “Modeling of ultrasonic wave propagation in teeth using PSpice: a comparison with finite element models,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 48, no. 4, pp. 1124–1131, 2001. View at Publisher · View at Google Scholar · View at Scopus
  26. N. N. Abbouda, G. L. Wojcikb, D. K. Vaughanb, J. Mould Jr., D. J. Powell, and L. Nikodym, “Finite element modeling for ultrasonic transducers,” in Medical Imaging: Ultrasonic Transducer Engineering, vol. 3341 of Proceedings of SPIE, San Diego, Calif, USA, February 1998. View at Publisher · View at Google Scholar
  27. R. Ramesh, C. D. Prasad, T. K. V. Kumar, L. A. Gavane, and R. M. R. Vishnubhatla, “Experimental and finite element modelling studies on single-layer and multi-layer 1–3 piezocomposite transducers,” Ultrasonics, vol. 44, no. 4, pp. 341–349, 2006. View at Publisher · View at Google Scholar · View at Scopus
  28. M. Mao, P. Verma, and W. J. Staszewski, “Finite element modeling for fault detection in medical ultrasonic transducers,” in Proceedings of the 5th European Workshop on Structural Health Monitoring, Sorrento, Italy, June-July 2010.
  29. J. Kocbach, Finite Element Modelling of Ultrasonic Piezoelectric Transducers, Influence of geometry and Material Parameters on Vibration, Response Functions and Radiated Field, University of Bergen Department of Physics, 2000.
  30. S. Moten, Finite Element Modeling of an Ultrasonic Transducer, Literature Survey, Eindhoven University of Technology, Eindhoven, The Netherlands, 2011.
  31. G. G. Yaralioglu, B. Baryram, A. Nikoozadeh, B. T. Khuri-Yakub, and E. L. Ginzton, “Finite element modeling of capacitive micromachined ultrasonic transducers,” in Proceedings of the Medical Imaging: Ultrasonic Imaging and Signal Processing, vol. 5750 of Proceedings of SPIE, San Diego, Calif, USA, February 2005. View at Publisher · View at Google Scholar
  32. T. P. Harrigan and E. E. Konofagou, “Estimation of material elastic moduli in elastography: a local method, and an investigation of Poisson's ratio sensitivity,” Journal of Biomechanics, vol. 37, no. 8, pp. 1215–1221, 2004. View at Publisher · View at Google Scholar · View at Scopus
  33. Q. Chen, S. I. Ringleb, T. Hulshizer, and K.-N. An, “Identification of the testing parameters in high frequency dynamic shear measurement on agarose gels,” Journal of Biomechanics, vol. 38, no. 4, pp. 959–963, 2005. View at Publisher · View at Google Scholar · View at Scopus
  34. K. Yoshikawa and G. Nakamura, “Model independent MRE data analysis,” Computational and Mathematical Methods in Medicine, vol. 2013, Article ID 912920, 11 pages, 2013. View at Publisher · View at Google Scholar · View at MathSciNet
  35. S. I. Ringleb, Q. Chen, D. S. Lake, A. Manduca, R. L. Ehman, and K.-N. An, “Quantitative shear wave magnetic resonance elastography: comparison to a dynamic shear material test,” Magnetic Resonance in Medicine, vol. 53, no. 5, pp. 1197–1201, 2005. View at Publisher · View at Google Scholar · View at Scopus
  36. Q. Chen, S. I. Ringleb, A. Manduca, R. L. Ehman, and K.-N. An, “A finite element model for analyzing shear wave propagation observed in magnetic resonance elastography,” Journal of Biomechanics, vol. 38, no. 11, pp. 2198–2203, 2005. View at Publisher · View at Google Scholar · View at Scopus
  37. Q. Chen, S. I. Ringleb, A. Manduca, R. L. Ehman, and K.-N. An, “Differential effects of pre-tension on shear wave propagation in elastic media with different boundary conditions as measured by magnetic resonance elastography and finite element modeling,” Journal of Biomechanics, vol. 39, no. 8, pp. 1428–1434, 2006. View at Publisher · View at Google Scholar · View at Scopus
  38. H. A. Naeeni and M. Haghpanahi, “Numerical study of shear wavelength observed in MRE experiments with FEM,” in World Congress on Medical Physics and Biomedical Engineering, September 7–12, 2009, Munich, Germany, vol. 25/4 of IFMBE Proceedings, pp. 700–703, Springer, Berlin, Germany, 2010. View at Publisher · View at Google Scholar
  39. H. A. Naeeni and M. Haghpanahi, “FE modeling of living human brain using multifrequency magnetic resonance elastograph,” Applied Mechanics and Materials, vol. 66–68, pp. 384–389, 2011. View at Publisher · View at Google Scholar
  40. Z. Hashemiyan, P. Packo, W. J. Staszewski, and T. Uhl, “Shear wave propagation modeling in magnetic resonance elastography using the local interaction simulation approach,” in Proceedings of the 11th World Congress on Computational Mechanics, 5th European Conference on Computational Mechanics and 6th European Conference on Computational Fluid Dynamics, Barcelona, Spain, 2014.
  41. J. L. Rose, Ultrasonic Waves in Solid Media, Cambridge University Press, Cambridge, UK, 1999.
  42. P. P. Delsanto, R. S. Schechter, and R. B. Mignogna, “Connection machine simulation of ultrasonic wave propagation in materials III: the three-dimensional case,” Wave Motion, vol. 26, no. 4, pp. 329–339, 1997. View at Publisher · View at Google Scholar · View at Scopus
  43. P. Paćko, T. Bielak, A. B. Spencer, W. J. Staszewski, T. Uhl, and K. Worden, “Lamb wave propagation modelling and simulation using parallel processing architecture and graphical cards,” Smart Materials and Structures, vol. 21, no. 7, Article ID 075001, pp. 1–13, 2012. View at Publisher · View at Google Scholar
  44. P. P. Delsanto, T. Whitcombe, H. H. Chaskelis, and R. B. Mignogna, “Connection machine simulation of ultrasonic wave propagation in materials. I: the one-dimensional case,” Wave Motion, vol. 16, no. 1, pp. 65–80, 1992. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  45. P. P. Delsanto, R. S. Schechter, H. H. Chaskelis, R. B. Mignogna, and R. Kline, “Connection machine simulation of ultrasonic wave propagation in materials. II: the two-dimensional case,” Wave Motion, vol. 20, no. 4, pp. 295–314, 1994. View at Publisher · View at Google Scholar · View at Scopus
  46. B. C. Lee and W. J. Staszewski, “Modelling of Lamb waves for damage detection in metallic structures. Part I. Wave propagation,” Smart Materials and Structures, vol. 12, no. 5, pp. 804–814, 2003. View at Publisher · View at Google Scholar · View at Scopus
  47. Z. Hashemiyan, M. Boeff, P. Packo et al., “Fault detection in medical ultrasonic transducers—comparison between finite element modeling and local interaction simulation approach,” Key Engineering Materials, vol. 588, pp. 157–165, 2012. View at Google Scholar
  48. P. Packo, T. Uhl, and W. J. Staszewski, “Generalized semi-analytical finite difference method for dispersion curves calculation and numerical dispersion analysis for Lamb waves,” The Journal of the Acoustical Society of America, vol. 136, no. 3, pp. 993–1002, 2014. View at Publisher · View at Google Scholar
  49. P. Packo, T. Bielak, A. B. Spencer et al., “Numerical simulations of elastic wave propagation using graphical processing units—comparative study of high-performance computing capabilities,” Computer Methods in Applied Mechanics and Engineering, vol. 290, pp. 98–126, 2015. View at Publisher · View at Google Scholar · View at MathSciNet
  50. J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 2nd edition, 2004. View at Publisher · View at Google Scholar · View at MathSciNet
  51. D. Iordache, P. P. Delsanto, and M. Scalerandi, “Pulse distortions in the FD simulation of elastic wave propagation,” Mathematical and Computer Modelling, vol. 25, no. 6, pp. 31–43, 1997. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus