Viscoelastic modulus reconstruction using time harmonic vibrations

This paper presents a new iterative reconstruction method to provide high-resolution images of shear modulus and viscosity via the internal measurement of displacement fields in tissues. To solve the inverse problem, we compute the Fr\'echet derivatives of the least-squares discrepancy functional with respect to the shear modulus and shear viscosity. The proposed iterative reconstruction method using this Fr\'echet derivative does not require any differentiation of the displacement data for the full isotropic linearly viscoelastic model, whereas the standard reconstruction methods require at least double differentiation. Because the minimization problem is ill-posed and highly nonlinear, this adjoint-based optimization method needs a very well-matched initial guess. We find a good initial guess. For a well-matched initial guess, numerical experiments show that the proposed method considerably improves the quality of the reconstructed viscoelastic images.


Introduction
Elastography [20] aims to provide a quantitative visualization of the mechanical properties of human tissues by using the relation between the wave propagation velocity and the mechanical properties of the tissues. During the last three decades, elastography led to significant improvements in the quantitative evaluation of tissue stiffness. The two major elastographic techniques are based on ultrasound [5,23,24,33,37] and on magnetic resonance imaging [17,18,34,20,21,22,25,26,30,32]. GE Healthcare has recently commercialized magnetic resonance elastography (MRE). Its main use is to assess mechanical changes in liver tissue. The mechanical properties of tissue include the shear modulus, shear viscosity, and compression modulus [14]. Quantification of the tissue shear modulus in vivo can provide evidence of the manifestation of tissue diseases. For centuries, palpation has been widely used to identify tissue abnormalities and estimate the mechanical properties of tissue. Therefore, it is surprising that the concept of remote palpation, which is the remote imaging of tissue stiffness, was first developed only in the late 1980s [9,12,27].
Although significant progress has been made in the development of shear modulus imaging technology, problems remain image quality relating to the enhancement of images of local tissue shear viscosity and shear modulus [10,13,15,19,29,30,31,35]. This paper focuses on the image reconstruction methods for tissue viscoelasticity imaging. To simplify the underlying inverse problem, the reconstruction of both the shear modulus and shear viscosity are considered under the assumption of isotropic elastic moduli.
This work considers the inverse problem of recovering the distribution of the shear modulus (µ) and shear viscosity (η µ ) from the internal measurement of the timeharmonic mechanical displacement field u produced by the application of an external time harmonic excitation at frequency ω/2π in the range 50 ∼ 200Hz through the surface of the subject. Modeling soft tissue as being linearly viscoelastic and nearly incompressible, the displacement u satisfies the elasticity equation where ρ denotes the density of the medium, ∇u t is the transpose of the matrix ∇u, λ is the compression modulus and η λ is the compression viscosity. The most widely used reconstruction method is the algebraic inversion method [17]: For any non-zero constant vector a, which requires the strong assumptions of ∇(µ + iωη µ ) ≈ 0 (local homogeneity) and (λ + iωη λ )∇ · u ≈ 0 (negligible pressure). The algebraic formula (2) ignores reflection effects of the propagating wave due to abrupt changes of µ + iωη µ , so that the method cannot measure any change of µ + iωη µ in the direction of a [13,28].
To deal with these fundamental drawbacks in the algebraic inversion method, the shear modulus decomposition algorithm based on Helmholtz-Hodge decomposition was developed in [13]. This is a much better performing method; however, it continues to neglect pressure by using (λ + iωη λ )∇ · u = 0, and is thus not realistic. In [29,30], the curl operator is applied to the elasticity equation (1) to eliminate the troublesome term (∇ × ∇((λ + iωη λ )∇ · u) = 0). The reconstruction method in [29,30] requires third-order derivatives of the noisy data u, making it very sensitive to noise in the data. A realistic model must take into account the non-vanishing pressure p [11,4], which can be defined roughly as p := lim λ→∞, ∇·u→0 (λ + iωη λ )∇ · u.
The shear viscoelasticity reconstruction method proposed in this paper is based on the full elasticity model. It does not require any derivative of u. The minimization of a misfit functional involving the discrepancy between the measured and fitted data is considered. The Fréchet derivatives of the functional with respect to µ and η µ are then computed by introducing an adjoint problem. This Fréchet derivatives based-iterative scheme requires a well-matched initial guess, because the minimization problem is highly nonlinear and may have multiple local minima. We find a well-matched initial guess that captures the edges of the image of the shear viscoelasticity.
The numerical results presented herein demonstrate the viability and efficiency of the proposed minimization method.

Viscoelastic model
Let an elastic subject occupy the smooth domain Ω ⊂ R d , d = 2, 3 with boundary ∂Ω. To evaluate the viscoelastic tissue properties, we create an internal time-harmonic displacement in the tissue by applying a time-harmonic excitation through the surface of the object. Under the assumptions of mechanical isotropy and incompressibility in the tissue, the induced time-harmonic displacement at angular frequency ω, denoted by u, is then governed by the full elasticity equation where ∇ s u = 1 2 (∇u + ∇u t ) is the strain tensor with ∇u t denoting the transpose of the matrix ∇u; ρ is the density of the medium; the complex quantity µ + iωη µ is the shear modulus, with µ indicating the storage modulus and η µ indicating the loss modulus reflecting the attenuation of a viscoelastic medium; λ and η λ are the compression modulus and compression viscosity, respectively. We assume that these heterogeneous parameters satisfy [14]: We define the interior domain Ω and the neighborhood E of the boundary, ∂Ω, as Ω := {x ∈ Ω|dist(x, ∂Ω) > with > 0}, E := Ω\Ω .
We assume that µ and η µ are known in the region E, and are denoted by µ 0 and η µ 0 , respectively. We denote by H s the standard Sobolev space of order s and by H s 0 the closure of C ∞ 0 , which is the set of C ∞ compactly supported functions, in the H s -norm.
Let us take Γ D ∪ Γ N = ∂Ω and Γ D ∩ Γ N = ∅. Boundary conditions on the displacement field u are imposed. Typically, we use an acoustic speaker system to generate harmonic vibration. If the acoustic speaker is placed on the portion Γ D of the boundary ∂Ω, then the boundary conditions for u can be expressed approximately by where n is the outward unit normal vector to the boundary.

Optimal control method
Define the misfit (or discrepancy) functional J(µ, η µ ) in terms of µ and η µ by the L 2 -norm in Ω of the difference between the numerical solution u[µ, η µ ] of the forward problem (30) and the measured displacement data u m = u m [µ * , η * µ ]: where µ * and η * µ are true distributions of shear elasticity and viscosity, respectively. The reconstruction of the unknowns µ and η µ can be obtained by minimizing the misfit functional J(µ, η µ ) with respect to µ and η µ .
In order to construct a minimizing sequence of J(µ, η µ ), we need to compute the Fréchet derivatives of J(µ, η µ ) with respect to µ and η µ . Assume that δ µ and δ ηµ are small perturbations of µ and η µ , respectively, by regarding δµ+iωδη µ µ+iωηµ ≈ 0. For notational simplicity, we denote u 0 := u[µ, η µ ], p 0 := the pressure corresponding to u 0 and p 0 + p 1 := the pressure corresponding to u[µ + δ µ , η µ + δ ηµ ]. Denoting the perturbation of displacement field by it follows from (30) that Let u 1 be the solution of the following problem Now we are ready to state two main theorems in this section which give the Fréchet derivatives of J(µ, η µ ) with respect to µ and η µ . Denote A : Furthermore, the Fréchet derivatives of J(µ, η µ ) with respect to µ and η µ are given by where v is the H 1 solution of the following adjoint problem [3,7]: The next theorem shows the differentiability of J(µ, η µ ).
In other words, if u 1 ∈ H 1 (Ω) is the weak solution to (8), as the perturbations δ µ , δ η → 0, we have the following formula: To prove the Fréchet differentiability Theorem 2.2 and the main Theorem 2.1, we need the following preliminary results.
Firstly, we state an interior estimate for the solution of the Stokes system whose proof basically follows from [6,8,16] by observing ∇ · ∇ s w = ∆w for w satisfying ∇ · w = 0. Lemma 2.3. For F ∈ L 2 (Ω) and (µ, η µ ) ∈ S, let w ∈ H 1 (Ω) be a weak solution of the following problem: Then, w ∈ H 2 (Ω) and where C is positive constant independent of F.
The following estimate for δu holds.
Proposition 2.4. The perturbation of displacement field δu ∈ H 1 (Ω) satisfies the following estimate: where C is positive constant independent of δ µ and δ ηµ .
Applying the interior estimate (12) to (13) and using Hölder's inequality and Sobolev embedding theorem [1,7], we arrive at This completes the proof. Now we are ready to prove Theorem 2.2.

Now we apply Proposition 2.4 to get
The proof is then completed. Now, it remains to identify the Fréchet derivatives of J(µ, η µ ). According to Theorem 2.2, the Fréchet derivatives ∂ ∂µ J(µ, η µ ) and ∂ ∂ηµ J(µ, η µ ) can be computed by expressing Ω u 1 (u 0 − u m )dx in terms of δ µ and δ ηµ . These are explained in the proof of Theorem 2.1.
Proof of Theorem 2.1. We use the adjoint solution v in (11) to get Using the vector identity ∇ · (qu 1 ) = ∇q · u 1 and divergence free conditions ( 0 = ∇ · δu = ∇ · u 1 = ∇ · v), the identity (16) can be rewritten as Since u 1 satisfies the equation (8), we have This proves the formula (9). The formula (10) can be obtained directly from Theorem 2.2 and the formula (9). This completes the proof.

Initial guess
Numerous simulations show that the reconstruction from an adjoint-based optimization method may converge to some local minimum that is very different from the true solution when the initial guess is far from the true solution. We observed that different initial guesses produce different reconstructions, and thus a good initial guess is necessary for accurate reconstruction using the iterative method (17).
We examine the optimization method using the initial guess obtained by the direct inversion method (2). Numerical simulations with this initial guess showed that serious reconstruction errors occur near the interfaces of different materials in the same domain; the direct inversion method cannot probe those interfaces. We found empirically that it is important to find an initial guess capturing the interfaces of different materials for the effective use of the optimization method.
To develop a method of finding such a good initial guess, we adopt the hybrid onestep method [15] which consider the following simplified model ignoring the pressure term: 2∇ · (µ + iωη µ )∇ s u + ρω 2 u = 0 in Ω, (18) where u is regarded as a good approximation of u[µ, η µ ]. To probe the discontinuity of (µ + iωη µ )∇ s u , we apply the Helmholtz decomposition where f and W are vector and matrix, respectively. The curl of matrix is defined in column-wise sense: where W j is the j-th column of matrix W for j = 1, 2, 3. Taking dot product of (19) with ∇ s u gives the following formula By taking the divergence to the equation (19), we have By taking the curl operation to the equation (19), we have Our proposed method for determining the initial guess is based on the modifying of hybrid one-step method. Using (21), an approximation of the vector potential f corresponding to the measurement u m can be computed by On the other hand, W can not be computed directly from u m since (22) contains unknown terms µ and η µ . Regarding µ + iωη µ in (22) as ∇ f :∇ sū m |∇ s um| 2 (see (20)), we can compute a rough approximation of W by solving Similarly, approximating µ + iωη µ by direct inversion formula (2), we can compute W by solving where a is any nonzero vector. Now, we use the formula (20) to get the initial guess of shear modulus by substituting f =f , W = (W 1 + W 2 )/2 and u = u m : In formula (26), the first term provides information in the wave propagation direction while the second term gives the information in the tangent direction of the wave propagation as shown in [15]. Note that if this initial guess is not satisfactory for the adjoint-based optimization problem, one can update the initial guess formula to obtain more accurate one by replacing (µ + iωη µ ) in (22) by (26). Numerical experiments demonstrates the possibility of probing the discontinuity of the shear modulus effectively. We emphasize that the initial guess plays an important role in Newton's iterative reconstruction algorithm based on the adjoint approach. By observing the adjoint problem (11), the load term u 0 − u m is related to the measured data and the initial guess in the first iteration step. If the initial guess ensure that ||u 0 − u m || is small in certain norm, the iteration scheme will converge and give good results. Otherwise, the initial guess makes ||u 0 − u m || far from 0 in certain norm, and the iteration scheme may not converge. This will be discussed in section 3.

Local reconstruction
In MRE, the time-harmonic displacement, u m , in the tissue is measured via phasecontrast-based MR imaging. Hence, the signal-to-noise ratio (SNR) of the data is related to that of the MR phase images, which varies from one region to another. For example, the SNR of data u m is very low in MR-defected regions, including the lungs, outer layers of bones, and some gas-filled organs. When the domain, Ω, contains such defected regions, the reconstructed image qualities may be seriously degraded by locally low SNR data in the defected regions. As a result, it would be desirable to exclude defected regions from Ω to prevent errors spreading in the image reconstruction. The proposed method is capable of a local reconstruction by restricting to a local domain of the interest. To be precise, let Ω loc be a subdomain of Ω in which u m has high SNR. Then, we consider the localized minimization problem with u loc [µ, η µ ] being the solution of in Ω loc , u = u m on ∂Ω loc .
As before, we need to compute the corresponding adjoint problem to get Fréchet derivative: There is no difference between the local reconstruction in Ω loc and the global reconstruction with Ω, except the boundary conditions. As in (17), the local reconstruction can be done by solving (28) and (29) with the initial guess (26). Local reconstruction requires that neither the boundary conditions need to be used on the whole domain, Ω, nor that the exact shape of Ω needs to be known. Numerical simulations verify the effectiveness of this local reconstruction, and will be discussed in section 3.

Numerical simulations
In this section, we perform several numerical experiments to illustrate the effectiveness of the shear viscoelasticity reconstruction algorithm proposed in the previous section.
To implement the reconstruction algorithm (17) proposed in section 2, we use the algorithm (20) in section 2.3 to initialize the iteration scheme. For numerical experiments, we set the two dimensional domain as Ω = [0, 10] × [0, 10] cm 2 with a boundary denoted by ∂Ω = Γ D ∪ Γ N ; see figure 3 (a). We apply the FEM method in Matlab (MathWorks In.) to solve the forward problem (30) as well as the adjoint problem (11) at each iteration step in the algorithm (17).
We set three different types of shear viscoelasticity distribution which are shown in the first column of figure 4 along with the true distribution of shear modulus and shear viscosity. The first and second rows are model 1, the third and fourth rows are model 2, and the fifth and sixth rows are model 3. For each model, the upper row shows elasticity while the lower row shows viscosity. Our numerical experiments are based on these three models. We generate two dimensional displacements u m = (u 1 , u 2 ) t by solving the problem (30) with frequency ω 2π =70Hz and density ρ = 1g · cm −2 . We apply the vibration to Γ D , and the other three sides boundaries are set to be traction free: For example, model 1 has the displacement fields shown in figure 3 where (b) and (c) are real parts of u 1 and u 2 , and (d) and (e) are imaginary parts of u 1 and u 2 , respectively. The next step is to implement our algorithm making use of these displacement fields with certain initial guesses of the distribution of viscoelasticity. We generate the initial guess by the direct inversion method (2) shown in the third column of figure 4 and the hybrid one-step method (20) shown in the fifth column of figure 4. From the generated initial guess, we can see that the reconstruction by the hybrid one-step method is much better than that of the direct inversion method in catching the inhomogeneous property of the medium. We have already explained the underlying mathematical reason for this phenomenon. We use the initial guesses from these two methods to initialize our proposed method, and the corresponding numerical results for each model are shown in the fourth column and last column of figure 4, respectively. For comparison, we also show the reconstruction with a homogeneous initial guess in each second column of figure 4.
The reconstruction results (see figure 4) show that the proposed method can reconstruct the viscoelasticity distribution with high accuracy (see (f) column) using a well-matched initial guess ( see (e) column). Otherwise, poor initial guesses (for example, the homogeneous initial guess and (c)), leads to unsatisfactory reconstructed images (see (b) and (d) columns). We also numerically evaluate the local reconstruction method proposed in section 2.4. We consider the rectangular domain, Ω, which is equally divided into four parts: top-left, top-right, bottom-left and bottom-right. It is assumed that the top-right part is contaminated by noise or defected data. For numerical simplicity, we add 3% white noise to the measured data in the top-right part. The reconstruction results in both the whole domain and the local domains are shown in figure 5 where (a) is the true distribution of shear viscoelasticity, (b) the initial guess with hybrid method, (c) the reconstruction in whole domain using proposed method, (d) the local reconstruction.

Conclusion
In this paper, we propose a reconstruction algorithm for shear elasticity and shear viscosity in a viscoelastic tissue. Our optimization-based approach involves introducing an adjoint problem to avoid taking any derivative of the measured time-harmonic internal data. The proposed initial guess formula is particularly suitable for imaging viscoelastic inclusions. The local convergence of the developed optimal control approach is an open problem. The recent stability results in [36] may be helpful in solving this difficult question. It would be also very interesting to generalize the proposed method for imaging anisotropic viscoelastic media. Another challenging problem is to recognize the disease state in tissue from multifrequency elastographic measurements. These important problems will be the subject of future work.