Estimation and validation of patient-specific liver elasticity distributions derived from 4DMR for radiotherapy purposes

Hepatocellular carcinoma (HCC) is one of the most common malignancies, and the third most common cause of cancer-related death worldwide (Jung et al 2013, Jun et al 2017). Surgical resection and liver transplantation are the primary treatment methodologies, however strict criteria limit the pool of eligible patients for both cases. HCC has a poor prognoses, with a 5-year survival rate of less than 12% due to a combination of late diagnosis and lack of efficient therapies for advanced stages (Affo et al 2017). Stereotactic body radiotherapy (SBRT) has been used to treat patients with HCC who are not eligible for other treatments (Feng et al 2017). SBRT uses advances in imaging and conformal radiotherapy to deliver ablative, high dose radiation in order to optimize local control. Radiation-induced liver disease (RILD) is a significant limiting factor in the use of SBRT because there are no effective treatments or predictors (Jung et al 2013). Most patients with HCC have pre-existing cirrhosis or hepatitis, which increases their risk of RILD (Kim et al 2015). Baseline liver function is thought to be the most important factor associated with risk of RILD (Jun et al 2017). Pre-treatment visualization of liver function in vivo is necessary in order to expand the use of SBRT for HCC.

The liver plays a role in metabolism, synthesis, secretion, immunity, and many other functions (Luna et al 2014). Liver disease, including cirrhosis, fibrosis and tumors, can disrupt the functions of the liver by altering the biomechanical properties of the tissue, most notably by changing the tissue stiffness (Sandrasegaran 2014, Li et al 2015). Clinically, elastography provides a measurement of liver stiffness and is a predictor for HCC (Pepin et al 2014). Magnetic resonance elastography (MRE) can noninvasively and quantitatively assess the elasticity characteristics of soft tissue (Li et al 2015). Liver stiffness measured by MRE has been shown to correlate well with histologic staging of fibrosis and differentiation of benign and malignant liver lesions (Venkatesh et al 2015). However, current MRE techniques require equipment that is not typically available within a radiotherapy setup.

Recently, an MR guided radiation therapy (MRgRT) system (MRIdian System™, ViewRay™, Cleveland, OH, USA) was introduced in the field of radiotherapy (RT) (Mutic and Dempsey 2014). The ViewRay system integrates a 0.35 T MRI with three Cobalt-60 sources, which allows for simultaneous MRI acquisition and treatment delivery for a variety of malignancies (Wojcieszynski et al 2016). The superior soft tissue contrast from Viewray's on-board MRI enables soft tissue based gated RT and on-line adaptive RT. Both techniques are especially helpful for treating thoracic and abdominal tumors by accounting for intra-fractional and inter-fractional motion. In this study, we present a novel elastography process for estimating liver tissue elasticity using 0.35 T 4DMRI to observe liver deformation in each respiratory cycle during patient simulation. The rest of this paper is organized as follows: section 2 introduces the methods employed for the deformable imaging registration, elasticity estimation, validation, and quantitative evaluation. Section 3 presents the qualitative and quantitative results of the deformable image registration, elastography, and validation processes. Section 4 presents a discussion of the results and highlights the areas for continued research, while section 5 concludes the paper.

Our elastography process focuses on estimating the effective Young's modulus for each voxel of liver tissue using 4DMR liver data. Figure 1 shows a flow chart summarizing the elasticity estimation. First, phase 1 and phase 8 datasets from the 4DMR liver images were registered using an optical flow deformable image registration (DIR) algorithm (Min et al 2014). Liver deformation vectors (DVFs) were obtained for every voxel of liver tissue. The biomechanical model was then assembled using segmented phase 1 liver geometry and a randomly initialized elasticity distribution. Using the ground-truth liver DVFs, the elasticity distribution was optimized. The inverse elasticity problem was formulated as a parameter-optimization problem with an objective to determine the elasticity parameter that would minimize the difference between the ground-truth deformation and the deformation computed by a biomechanical model. The biomechanical model and inverse elasticity estimation was implemented on a GPU cluster, which allowed the elasticity estimation for each patient dataset to converge in around 2 h. Spatial elasticity and displacement error distributions were then obtained and validated.

Zoom In Zoom Out Reset image size

In this section, we first briefly present the 4DMR acquisition process and the DIR technique used to obtain the ground-truth DVF. Next, we present the constitutive model and convergence criteria used for the elasticity estimation process. We then describe a technique that was used to validate the estimated elasticity. We conclude the section with a discussion of the metrics used to evaluate the quantitative accuracy of the elasticity estimation and validation results.

2.1. 4DMRI acquisition

2.2. Deformable image registration

2.3. Elasticity estimation

The inverse elasticity problem was formulated as a parameter-optimization problem with an objective to determine the elasticity parameter that would minimize the difference between ground-truth DVFs and those computed by a constitutive biomechanical model. Solving the inverse elasticity problem was carried out in two steps: (1) estimating the DVF for every voxel of liver tissue for a given elasticity distribution and boundary constraint using the biomechanical model, and (2) optimizing the elasticity distribution to best reproduce the ground-truth DVF. These steps will be further explained in the following sections.

2.3.1. Biomechanical model

The first step requires a forward biomechanical constitutive model. The model presented in this section is focused on computing the liver tissue deformation for a give boundary constraint arising from the ground-truth DVF. This model has been previously described for the head and neck, breast, and lung anatomies in (Neylon et al 2015, Hasse et al 2016, Hasse et al 2017a, Hasse et al 2017b) but will be summarized briefly here. This approach is appropriate for our current study since its applicability has been systematically investigated for breathing motion, which is the main factor in diaphragm-induced liver deformations (Wu et al 2008).

The biomechanical liver model was assembled from the segmented 4DMR phase 8 dataset. The model geometry was represented by finite element nodes corresponding to the center of each voxel of the phase 8 image. Liver tissue was viewed as a series of finite mass elements coupled with linear elastic connections to adjacent mass elements in deformation space to ensure a physically realistic deformation. Our approach uses linear elastic constitutive laws that were implemented because of the quasi-static tissue response of liver tissue to small diaphragm-induced deformations (Nava et al 2008).

The model is actuated by enacting changes in boundary constraints, which causes corrective forces to be applied to each mass element. The boundary constraints were computed from the liver boundary DVF. The corresponding corrective forces on each mass element were calculated by summing elastic, shear, and dashpot damping forces for each connected elements. The elastic force $({\overrightarrow{f}}_{E,ab})$ acting between two mass elements a and b is described by:

$$\begin{eqnarray}&&{\overrightarrow{f}}_{E,ab}=\displaystyle \sum _{b}\left({E}_{ab}* \displaystyle \frac{{\rm{\Delta }}{L}_{ab}}{{L}_{ab}}\right),\end{eqnarray} \tag{ 1 }$$

where E_ab is the effective Young's modulus acting between the two elements a and b, and L_ab is the rest length distance, and ΔL_ab is the change in length between the two elements. Since liver tissue is modeled as a linearly elastic material, (1) is derived from the inverse relationship between Young's modulus and displacement; for a constant force, if elastic modulus increases then displacement must decrease and vice versa (Doyley et al 2000).

The shear spring force $({\overrightarrow{f}}_{S,ab})$ on element a due to element b is described by

$$\begin{eqnarray}&&{\overrightarrow{f}}_{S,ab}=\displaystyle \sum _{b}\left({S}_{ab}* \displaystyle \frac{{L}_{ab}-{\rm{\Delta }}{L}_{ab}}{{L}_{ab}}\right),\end{eqnarray} \tag{ 2 }$$

where S_ab is the shear moduli between elements a and b, set to 4 kPa for each voxel in our estimations, an approximation based on the relationship between shear modulus, Young's modulus, and Poisson's ratio and similar to values seen in (Liu et al 2010, Suki et al 2011). The dashpot damping force is calculated from the relative velocities $\overrightarrow{v}$ of elements a and b, and a local damping factor μ_ab

$$\begin{eqnarray}&&{\overrightarrow{f}}_{v,ab}=\displaystyle \sum _{b}({\mu }_{ab}* (\overrightarrow{{v}_{b}}-\overrightarrow{{v}_{a}})).\end{eqnarray} \tag{ 3 }$$

The new positions and velocities of each mass element $({\overrightarrow{x}}_{a}^{n+1},{\overrightarrow{v}}_{a}^{n+1})$ were then updated from the previous values $({\overrightarrow{x}}_{a}^{n},{\overrightarrow{v}}_{a}^{n})$ using Implicit Euler integration and the total internal corrective force $\overrightarrow{{f}_{a}}.$

$$\begin{eqnarray}&&{\overrightarrow{v}}_{a}^{n+1}={\overrightarrow{v}}_{a}^{n}+\left(\displaystyle \frac{{\overrightarrow{f}}_{a}}{{m}_{a}}+\overrightarrow{g}\right)\delta ,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\overrightarrow{x}}_{a}^{n+1}={\overrightarrow{x}}_{a}^{n}+{\overrightarrow{v}}_{a}^{n+1}\delta .\end{eqnarray} \tag{ 5 }$$

In (4) and (5), δ is the time step between iterations, m_a is the mass of each mass element a, and $\overrightarrow{g}$ is the acceleration due to gravity. The distance between the new and previous positions of each mass element was taken to be the displacement for that mass element, henceforth denoted as d_a.

2.3.2. Optimization process

The iterative scheme for estimating the elasticity distribution is based on the inverse relationship between Young's modulus and displacement. Given an initial elasticity distribution, the constitutive model described above computed the displacement of every mass element. These displacements were compared with the ground-truth DVF values, and the elasticity distribution was updated until the differences was minimal. New displacements were generated for each iteration by updating the effective elasticity distribution as described below.

Initial (E_a), minimum (E_min,a), and maximum (E_max,a) Young's modulus values were initialized for each voxel a in the biomechanical model. For liver elastography purposes, the Young's modulus range was set to be from 1 kPa to 25 kPa, an extended approximation of values from the literature (Nitta et al 2009). The displacement differential was calculated according to the following relationship:

$$\begin{eqnarray}&&{\rm{\Delta }}{d}_{a}=\parallel {d}_{0,a}-{d}_{a}\parallel ,\end{eqnarray} \tag{ 6 }$$

where d_0,a and d_a refer to the ground-truth and model generated displacements of voxel a, respectively. The displacement differential was then used to binaurally update the E_min,a and E_max,a:

$$\begin{eqnarray}&&if\,{\rm{\Delta }}{d}_{a}\lt 0:\,\left\{\begin{array}{l}{E}_{\min ,{a}}={E}_{a}\\ {E}_{\max ,{a}}={E}_{\max ,{a}}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&if\,{\rm{\Delta }}{d}_{a}\gt 0:\,\left\{\begin{array}{l}{E}_{\min ,{a}}={E}_{\min ,{a}}\\ {E}_{\max ,{a}}={E}_{a}\end{array}\right..\end{eqnarray} \tag{ 8 }$$

A new E_a', which became the initial value for the next iteration, was then generated at each iteration according to the modified Gauss-Newton optimization scheme (Fu et al 2000) with a step size:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{a}=({E}_{\max ,{a}}-{E}_{\min ,{a}}),\end{eqnarray} \tag{ 9 }$$

giving us:

$$\begin{eqnarray}&&{E}_{a}^{{\prime} }={E}_{a}+{\rm{\Delta }}{E}_{a}.\end{eqnarray} \tag{ 10 }$$

After the elasticity was updated, the biomechanical model was then reset to the initial un-deformed rest position, and the boundary constraints were re-applied. The reconstruction process was repeated until a suitable stopping (convergence) criteria was reached.

2.4. Quantitative evaluation and validation

The convergence metrics utilized in this study are critical to interpreting the results. At each iteration of the elasticity estimation, the resultant model DVF was compared to the ground-truth DVF. The number of voxels that converged with a certain tolerance was quantified as follows:

$$\begin{eqnarray}{c}_{a}=\left\{\begin{array}{lll}1 & : & \parallel | {d}_{a}| -| {d}_{a,0}| \parallel \lt {\epsilon }\\ 0 & : & otherwise\end{array}\right.\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{p}_{{\epsilon }}=\displaystyle \frac{\displaystyle {\sum }_{1}^{n}{c}_{a}}{n}.\end{eqnarray} \tag{ 12 }$$

In (11), represents the L2-norm of the DVF for each voxel a and c_a represents the error for each voxel. The percent accuracy, ${p}_{{\epsilon }}$ was then tabulated by summing the c_a over all voxels. For elastography purposes, the iterative process was stopped when ${p}_{{\epsilon }}\geqslant 0.95$ or after 100 iterations. The convergence criteria was formulated to be 1 mm, the current standard for DIR-based targeting errors for liver SBRT dose planning purposes (Velec et al 2011).

An additional error metric was necessary in order to ensure the precision of the estimated elasticity. Specifically, we needed to ensure that the estimated elasticity and displacement distributions could be used to warp the phase 1 liver dataset to closely represent the phase 8 image. Furthermore, the accuracy of the ground-truth data needed to be more fully quantified. For this purpose, normalized cross-correlation was used as an imaging metric to quantify accuracy between both the target and warped images along with the model-deformed and ground-truth deformed images. Normalized cross-correlation was used because it has been found to be an accurate similarity measurement for the registration of soft tissue, and is unaffected by intensity changes (Penney et al 1998).

A key part of our focus is to validate the estimated elasticity using clinical datasets. However, patient-specific elasticity distributions have not been previously documented in the literature. For validation purposes, deformation vectors from the phase 4 datasets (representing mid-expiration breathing phase) that were not used for the elasticity estimation were employed. The accuracy of the elasticity distribution derived from phase 1 and phase 8 datasets was quantified by its' ability to also represent the phase 4 liver data. To further ensure the quantitatively evaluate the precision of the DIR and the estimated elasticity, normalized cross-correlation (NCC) was used as an image similarity metric between two sets of images for both the estimation and validation process: first between the target and warp images used in the deformable image registration, and second between the target image and model-deformed images.

In this section, we first present the qualitative DIR results, followed by qualitative elasticity estimation and validation results. Finally, we present the quantitative results for the respective sections.

3.1. Deformable image registration

The accuracy of the DIR results is necessary to ensure the validity of the ground-truth deformation data. Figure 2 illustrates an example of the DIR for a 2D-slice of 4DMR data. Figure 2(a) shows the source data (phase 8 dataset) in red overlain with the target (phase 1) data in green. Mismatches can be seen between the source and target, especially near the diaphragm and liver boundaries. Figure 2(b) shows the target data is red overlain with the warped source data in green. It can be seen that the registration accounts for the liver deformation well, with little mismatch seen between the target and warp images. The extent of the registration accuracy will be further quantified in section 3.3.

Zoom In Zoom Out Reset image size

3.2. Elasticity estimation.

Figure 3 shows the displacement accuracy corresponding with a resultant elasticity map for a 2D-slice of patient liver. Figure 3(a) shows the source image of a segmented liver. Figure 3(b) shows the ground-truth displacement derived from DIR in mm, while figure 3(c) shows the optimized model displacement, also in mm. Figure 3(d) shows the displacement error between 3(b) and 3(c), with error greater than 1 mm highlighted in blue. Finally, figure 3(e) shows the optimized elasticity distribution. The low error supports the validity of the error estimation. Furthermore, regions of high elasticity correspond to apparent liver lesions seen in figure 3(a).

Zoom In Zoom Out Reset image size

3.3. Quantitative evaluation and validation

Table 1 indicates the average elasticity, displacement, and accuracy for both the estimation and validation for each patient dataset. Overall, 95% of voxels converged within 1 mm of ground-truth. The average elasticity ranged from 2.69 kPa to 6.42 kPa. The maximum deformation ranged from 3.73 mm to 9.04 mm. Variations in the liver physiological factors, such as range of motion, did not affect the accuracy, indicating the applicability of the elastography process to a wide range of patients.

Table 1.  Elasticity estimation and validation results for 11 patient datasets.

	Estimation			Validation
Patient	Elasticity (kPa)	Maximum displacement (mm)	Error < 1 mm	Maximum displacement (mm)	Error < 1 mm (%)
1	5.88	9.04	71.34	7.88	53.10
2	2.69	5.10	98.88	5.08	97.94
3	5.91	4.96	97.61	3.05	86.94
4	6.16	6.29	94.64	3.76	91.42
5	4.90	6.30	98.97	5.83	85.61
6	6.26	8.60	99.13	6.56	88.22
7	6.42	6.61	96.17	5.71	95.07
8	5.80	7.64	96.35	5.73	69.92
9	4.44	3.73	100	2.79	100
10	5.28	3.99	99.98	2.59	99.584
11	5.06	6.34	96.61	6.20	76.67
AVG	5.34	6.23	95.43	5.02	85.86

As no ground-truth elasticity distributions are available for absolute elasticity accuracy calculations, validation was performed using phase 4 (representing mid-expiration breathing phase) 4DMR liver data. To illustrate the validation process and results, we present figure 4. Figure 4(a) shows the elasticity estimated from the ground-truth displacement vectors. Figure 4(b) shows the validation displacement acquired from deformable image registration between phase 1 and phase 4 images. Figure 4(c) shows the resultant model displacement in mm after the liver biomechanical model is deformed with the elasticity seen in figure 4(a). Figure 4(d) shows the displacement error between (b) and (c), with error greater than 0.3 mm highlighted in blue. A majority of the voxels converged with high accuracy, illustrating the estimated elasticity's validity in representing different 4DMR phases. The validation results are quantitated in table 1; maximum deformation ranges from 2.59 to 7.88 mm, and average accuracy is 86%.

Zoom In Zoom Out Reset image size

Finally, image similarity results are presented in table 2 below. On average, the deformable image registration had an accuracy of around 0.97 for both the estimation and validation processes. The resultant model-deformation represented the ground-truth deformation with a similarity of 0.97 for the estimation process and 0.90 for the validation process.

Table 2.  Image similarity metric normalized cross-correlation between target (phase 1) and warp (deformed phase 8) data and optimized model deformation and ground-truth deformation for both estimation and validation results.

	Estimation		Validation
	Target → Warp	Model DVF → Ground-truth DVF	Target → Warp	Model DVF → Ground-truth DVF
1	0.941	0.978	0.952	0.866
2	0.966	0.861	0.967	0.779
3	0.982	0.972	0.977	0.931
4	0.979	0.944	0.977	0.795
5	0.974	0.993	0.973	0.947
6	0.972	0.982	0.977	0.949
7	0.974	0.953	0.972	0.896
8	0.973	0.994	0.984	0.912
9	0.963	0.993	0.965	0.952
10	0.954	0.995	0.961	0.916
11	0.966	0.990	0.973	0.971
Average	0.968	0.969	0.971	0.901

In this paper, we presented the results of a liver elastography process performed on 11 4DMR datasets. A physics-based biomechanical model was used to solve the inverse elasticity problem. Liver DVFs from the registration of phase 1 and phase 8 diaphragm positions were obtained using an in-house optical flow DIR algorithm. Liver boundary displacements were employed as boundary constraints, while the inner liver tissue voxels were allowed to deform according to linear elastic material properties.

On average, 95% of voxels for 11 patients converged within 1.0 mm of ground-truth deformation. Maximum deformation ranged from 3.73 to 9.04 mm for the estimation cohort. The average elasticity ranged from 2.69 to 6.42 kPa. The average values found here correspond well with those found in the current work in the field (Singh et al 2015, Venkatesh et al 2015, Zeng et al 2017). In addition, an image similarity metric showed high similarity between the phase 8 and warped phase 1 registration and experimental results, with values of 0.97 and 0.97 respectively. Overall, these results suggest that liver elastography can be performed using ViewRay 4DMR datasets for a wide range of patients.

The potential of 4DMR liver elastography to be used in the clinic requires extensive validation. Phase 4 datasets were obtained as a validation cohort so that the elasticity results could be validated in a clinically-relevant manner. The maximum deformation ranged from 2.59 to 7.88 mm for the validation cohort, and 86% of voxels converging within 1 mm of clinical ground-truth deformation. An image similarity metric again showed high similarity between phase 4 and warped phase 1 registration and validation model results, with values of 0.97 and 0.90 respectively. Future work will investigate obtaining MRE data for patients so the elasticity distributions can be more explicitly and quantitatively validated.

We foresee two main limitations to the methodology presented here. First, the linear elastic constitutive model employed in the biomechanical model might be overly simplistic for describing the diaphragm-induced motion of diseased liver tissue. Future studies will focus on expanding to a non-linear constitutive model. Secondly, since the ViewRay implementation at UCLA is relatively new, we only had access to a small number of patients. Future work will focus on collecting more patients for further studying the elastography process along with a multi-modal validation.

The 95% convergence of the elastography technique presented here indicates the potential for integrating liver elastography directly within the radiotherapy workflow. Elastography of the liver anatomy using (Viewray) 0.35 T MR Scanners has not yet been developed and validated. The elastography process described in this manuscript can be implemented within the radiotherapy treatment setup without the need for additional imaging modalities. Elastography information acquired within the radiotherapy setup can lead to improved tissue sparing radiotherapy treatment plans and more precise monitoring of treatment response.

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1144087, the US Department of Defense Virtual Tissue Consortium, and the UCLA Department of Radiation Oncology.