Journal of Healthcare Engineering

Volume 2019, Article ID 2374645, 12 pages

https://doi.org/10.1155/2019/2374645

## An Iterative Method for Estimating Nonlinear Elastic Constants of Tumor in Soft Tissue from Approximate Displacement Measurements

^{1}Department of Biomedical Engineering, Amirkabir University of Technology (AUT), Tehran 15875-4413, Iran^{2}Faculty of Medical Sciences and Technologies, Science and Research Branch, Islamic Azad University, Tehran 1477893855, Iran

Correspondence should be addressed to Ali Fallah; ri.ca.tua@hallafa

Received 24 March 2018; Revised 24 June 2018; Accepted 12 July 2018; Published 6 January 2019

Academic Editor: Zahid Akhtar

Copyright © 2019 Maryam Mehdizadeh Dastjerdi 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.

#### Abstract

*Objectives*. Various elastography techniques have been proffered based on linear or nonlinear constitutive models with the aim of detecting and classifying pathologies in soft tissues accurately and noninvasively. Biological soft tissues demonstrate behaviors which conform to nonlinear constitutive models, in particular the hyperelastic ones. In this paper, we represent the results of our steps towards implementing ultrasound elastography to extract hyperelastic constants of a tumor inside soft tissue. *Methods*. Hyperelastic parameters of the unknown tissue have been estimated by applying the iterative method founded on the relation between stress, strain, and the parameters of a hyperelastic model after (a) simulating the medium’s response to a sinusoidal load and extracting the tissue displacement fields in some instants and (b) estimating the tissue displacement fields from the recorded/simulated ultrasound radio frequency signals and images using the cross correlation-based technique. *Results*. Our results indicate that hyperelastic parameters of an unidentified tissue could be precisely estimated even in the conditions where there is no prior knowledge of the tissue, or the displacement fields have been approximately calculated using the data recorded by a clinical ultrasound system. *Conclusions*. The accurate estimation of nonlinear elastic constants yields to the correct cognizance of pathologies in soft tissues.

#### 1. Introduction

According to the World Health Organization report, cancer is one of the principal morbidity and mortality agents throughout the world, with approximately 8.8 million deaths (nearly 1 in 6 deaths) in 2015 and 70% increase in the number of new cases over the next two decades. The cancer statistics imply the requisite to extend medical scrutiny to improve cancer prevention, early correct diagnosis, meticulous screening, and effective treatment and reduce the invasiveness and costs of applied techniques.

Since the first introduction of ultrasound (US) imaging in clinical practice in the 1970s, ultrasonography and other US modalities, for example, Doppler imaging and state-of-the-art elastography imaging methods, which provide the information related to the tissue acoustic impedance, vascular flow, and tissue mechanical characteristics or variables such as its stiffness or strain, respectively, have been extensively utilized for medical diagnoses [1]. The US imaging is recognized a noninvasive, safe, easy-to-use, low-cost, and widely accessible imaging modality for visualizing in vivo tissues. Elastography approach has currently been regarded a promising alternative to invasive medical procedures, for example, the biopsy, to characterize tissue abnormalities.

The wide variety of strategies that are being employed to quantify and image mechanical properties of biological tissues are recognized as elastography or elasticity imaging techniques with reference to their similar premise [2]:(1)The in vivo tissue is being deformed by a specified external or internal load or motion.(2)The response of tissue is being recorded by the use of a standard clinical imaging system, such as the US or magnetic resonance imaging (MRI) system.(3)The mechanical characteristics of tissue are estimated through the assessment of tissue displacement fields.

The alterations in the microstructure of tissue as a consequence of pathophysiological phenomena would change the mechanical properties of tissue; for instance, the increase in the stromal density of cancerous tissue would cause the increment in its Young’s modulus [3, 4]. The outcomes of numerous experimental research studies carried out by Krouskop et al. [5], Samani et al. [6, 7], Lyshchik et al. [8], Soza et al. [9], Hoyt et al. [10], Schiavone et al. [11], O’Hagan et al. [12], and Moran et al. [13], to mention but a few, have confirmed the relation between tissue structures and macroscopic mechanical features which are being evaluated, quantified, and/or imaged by employing the palpation, elastography, digital rectal examination, and such like methods.

The precise determination of mechanical characteristics of understudy tissue by the use of an elastography technique would undoubtedly necessitate realistically appraising or modeling the tissue manners, specifically its nonlinear response to the stimulation. Hyperelasticity theory is one of the constitutive theories that have been manipulated to model the nonlinear constitutive demeanor demonstrated by biological soft tissues. A variety of hyperelastic models, for instance the well-known Neo-Hookean, Mooney–Rivlin, Yeoh, and Polynomial models, have been recommended for this purpose [14–17]. In comparison with the studies involving the linear elasticity imaging, the number of research studies conducted to image the nonlinear features of tissues is limited [2].

With the aim of diagnosing a tumor inside the understudy tissue correctly, we have utilized an iterative method, as explicated in the next section, to accurately estimate the Mooney–Rivlin hyperelastic parameters of the tumor inside the tissue. The displacement fields inside the tumor have been analyzed to extract its hyperelastic parameters. An iterative technique has been employed since it has been assumed that no initial knowledge of the tumor was accessible except the displacement fields inside the tumor. The displacement fields inside the tumor have been extracted from the simulated/recorded radio frequency (RF) signals using the cross correlation-based method. The response of an abnormal tissue to a sinusoidal load (with low frequency to negate the inertia) has been simulated by applying the finite element (FEM) software package ABAQUS. The RF signals have been simulated by the use of the Field II US Simulation Program. In brief, in this paper, we scrutinize the diagnosis of tumor through its hyperelastic parameters in the conditions where there is no primary perception of the tumor and the displacement fields inside the tumor are estimated imprecisely.

#### 2. Materials and Methods

##### 2.1. Hyperelasticity Theory

Constitutive theories take the improvement of mathematical models, also known as constitutive equations, into consideration in order to provide the possibility to minutely describe the behavioral characteristics of materials. Constitutive theories in continuum mechanics deal with formulating material models that are [18, 19](a)on the basis of some mechanical universal principles(b)in accordance with experimental observationsThoughtful consideration of soft tissue’s model would result in the realistic prediction of its behavior such that it could be verified by experimental observations. Nonlinear constitutive manners that have been observed from soft tissues in numerous in vivo and ex vivo experimental research studies could be modeled by the use of hyperelastic models [15, 16, 20].

The hyperelastic constitutive laws deal with modeling materials with nonlinear elastic behaviors in reaction to large strains. The nonlinearities that are the consequences of (a) the material behavior and (b) the significant change in the shape of material are both regarded in the constitutive theory of hyperelastic materials. Hyperelastic materials are generally described by specific forms of strain energy density (stored energy) functions. While characterizing the homogeneous material’s absorbed energy due to its deformation, the strain energy function, , is defined as a function of deformation gradient, [21, 22]

It is considered that B_{r} and B represent, respectively, the reference or undeformed configuration, which refers to the situations where no load is exerted to the material and the deformed configuration, which is relevant to the situations where the material is under load and therefore it may alter with time, *t*. In addition, it is assumed that **X** and **x**, respectively, correspond to the position vectors of a material point in the reference and deformed configurations, B_{r} and B. The time-dependent deformation of material, that is, the motion of material point, from B_{r} to B could be described by the function **χ**, which, for each *t*, (a) is an invertible function and (b) satisfies proper regularity conditions as follows [23]:

The deformation gradient tensor, **F**, is defined aswith Cartesian componentswhere Grad, *x*_{i}, and *X*_{α} refer to the gradient operator in the configuration B_{r} and components of **x** and **X**, respectively, while the general convention,is satisfied. Due to the local invertibility of deformation, **F** should be nonsingular. The unique polar decomposition of **F** is defined aswhere the tensor **R** is an appropriate orthogonal tensor and the tensors **U** and **V** are symmetric positive-definite tensors known as the right and left stretch tensors, respectively. Equation (7) represents the spectral decompositions of tensors **U** and **V**:where each *λ*_{i} refers to one of the principal stretches, and are the unit eigenvectors of and known as the Lagrangian and Eulerian principal axes, and is the sign of tensor product [23, 24].

The tensor function is regarded as the function of material response in the configuration B_{r} with respect to the nominal stress, , that is, the transpose of the first Piola-Kirchhoff stress. The following equation for the nominal stress, ,is validated for an unconstrained homogeneous hyperelastic material. While the material is incompressible, the arbitrary hydrostatic pressure, *p*, which is the Lagrange multiplier associated with the material incompressibility, modifies the relation for the nominal stress, , as

With regard to the relation between the nominal stress tensor, **S**, and Cauchy stress tensor, ,the Cauchy stress tensor, , could be calculated through (11) in which the symmetric tensor function **G** denotes the function of material response in the configuration B_{r} associated with the Cauchy stress tensor, ,

With respect to the arbitrary hydrostatic pressure, *p*, defined previously, for an incompressible material, the relation for the Cauchy stress tensor, , modifies as [23–25]

First, second, and third invariants of , known as the strain invariants of deformation, which make provision for mapping the area and volume between the deformed configuration, B, and reference configuration, B_{r}, are computed throughfor an unconstrained isotropic elastic material. The left Cauchy-Green deformation tensor, **B**, and its principal invariants are calculated as follows:

For incompressible materials, a slightly different set of principal invariants of **B**, as represented, is generally employed:

The right Cauchy–Green deformation tensor, **C**, and its principal invariants are computed similarly.

The Cauchy stress tensor for an unconstrained isotropic elastic material is computed in terms of strain invariants, *I*_{1}, *I*_{2}, and *I*_{3}, as follows:as the result of the absence of (because *I*_{3} = 1) and the presence of *p* for incompressible materials, the Cauchy stress tensor changes tofor the forenamed materials [21–24], which simplifies to

A variety of stored energy functions have been introduced in the literature that could be employed to model the nonlinear elastic behavior of soft tissues precisely [26, 27], from the popular long-standing Neo-Hookean model (originated by Treloar in 1943 [28]) and Mooney–Rivlin model (proposed by Rivlin et al. in 1951 [29, 30]) to the state-of-the-art models, for instance, the ones introduced by Limbert in 2011 [31], Nolan et al. in 2014 [32], and Shearer in 2015 (for modeling ligaments and tendons) [33], although this is hitherto an active field of study in the material and biomedical sciences.

Being pertinent to model the behavior of a wide range of materials, for instance the soft tissues and polymers, the Mooney–Rivlin model is one of the conventional hyperelastic models in the literature [34–39]. Two historical attributes have made this model a distinguished one: (1) it is one of the primarily introduced hyperelastic models; (2) it could meticulously predict the nonlinear demeanor observed from some materials, specifically the isotropic rubber-like materials. The most general version of this model has been defined based on (the linear combination of) the first and second strain invariants of deformation, and . The Mooney–Rivlin strain energy function is expressed aswhere and are the material hyperelastic constants, is the constant related to the material volumetric response (i.e., the material bulk modulus), and *J* is the determinant of deformation gradient tensor, **F**. In addition, regarding the initial shear modulus of material, *μ*_{0}, the relationlinks the two hyperelastic parameters [26, 36]. For incompressible materials, the Mooney–Rivlin strain energy function simplifies to (since *J* = 1)

##### 2.2. Soft Tissue Simulation

To estimate nonlinear elastic parameters of an unidentified tumor using the proposed technique, we have simulated a simplified 3D breast tissue geometry utilizing the FEM software ABAQUS (Dassault Systèmes Simulia Corp., Johnston, RI, USA). The breast tissue is comprised of three partial tissues, fat, fibroglandular, and tumor, located consecutively from outside to inside. We have applied the Mooney–Rivlin hyperelastic model to evaluate the deformation of simulated breast tissue induced by the external excitation.

For estimating hyperelastic parameters of the tumor, we applied a sinusoidal load with frequency of 0.1 Hz (the very low frequency to ignore the inertia effects) to the simulated tissue and registered its response to extract the displacement fields inside the tumor. It is feasible to estimate displacement quantities inside the in vivo tissue by the use of some conventional medical imaging systems such as the US imaging or MRI system. We discuss the methods of calculating the displacements inside the tissue from RF signals or images recorded by the US imaging system in Section 2.4. In order to provide the accessibility to displacement values at some sequential moments, the US images or RF signals should be continuously saved for a period of time.

##### 2.3. Simulation of US RF Signals and Images

When a tissue is inspected with the US imaging system, the tissue is scanned with respect to the probe; in other words, when the tissue is compressed with the probe, the presented features of the tissue in the image seemingly move upward, although in point of fact they would move downward. In furtherance of simulating the postcompression RF signals and B-mode images in the probe coordinate system, instead of the phantom’s surface in contact with the probe, the surface in front of the probe was assigned to be moving [40–42].

The RF signals and B-mode images of the simulated phantom while responding to the sinusoidal load have been simulated using the Field II US Simulation Program (A MATLAB® toolbox for US field simulation) [43]. The nodal displacement measurements of the phantom in response to the sinusoidal load have been applied to simulate the postdeformation RF signals and B-mode images in varied deformation states. With the view to simulating the RF signals and B-mode image correlated with a particular deformation state, the correspondent nodal displacement values achieved by the finite deformation analysis have been linearly interpolated to compute the positions of scatterers corresponding to the specified deformation state.

##### 2.4. Estimation of Displacement Field

In addition to the US elastography imaging techniques, motion tracking algorithms have been applied in various US-based methods, for example, blood flow imaging, thermal strain imaging, phase-aberration correction, strain compounding, and temperature imaging, to name a few. The prominence of clinical applications of motion tracking methods has contributed to the significant accrual in the number of relevant investigations and the proposal of a multitude of techniques including phase-domain tactics, time-domain (1D) or space-domain (2D) procedures, and spline-based methods [44].

Among the propounded techniques, the cross-correlation algorithm is known as the gold standard of motion estimation. In this technique, the displacement quantities are computed through searching the locations of the maximums of cross-correlation values between corresponding axial-lateral grids of small windows (with high overlap) in the pre- and postdeformation frames recorded from the medium. The shifts between the pre- and postdeformation windows quantify the displacement field in the medium [44, 45].

The cross-correlation algorithm with guided search has been applied to initially estimate the displacement quantities in the tumor at the selected step times from the simulated pre- and postdeformation RF signals and B-mode images, and afterwards evaluate the errors of displacement estimates. The calculated displacement values of proximate regions (i.e., previous samples and lines) have been exploited to reduce the search area and therefore significantly decrease the computational expense of the cross-correlation algorithm. By the use of simulated RF signals, although we devoted more time to calculate the displacement values (because of their higher sampling rate) compared with the US images, we achieved more accurate displacement estimates with higher spatial resolution.

##### 2.5. Estimation of Hyperelastic Parameters

Regarding the aforementioned explanation related to the finite strain theory (also known as large deformation theory), when a uniaxial stress, **σ**, is applied to the material, the deformation gradient tensor, **F**, could be computed aswhere *λ*_{1}, *λ*_{2}, and *λ*_{3} are the principal stretches with respect to the set of coordinate axes (corresponding with *x*_{i} = *λ*_{i}*X*_{i}, *i* = 1,2,3). The principal invariants, *I*_{1}, *I*_{2}, and *I*_{3}, could be stated further in terms of principal stretches as follows:

If **λ** symbolizes the stretch parallel to the stress applied to the medium in line with the first coordinate axis (that is equal to the *λ*_{1} set), two assumptions, (1) the equality of deformations in the two other coordinate axes and (2) the incompressibility of the medium (*I*_{3} = 1), result in simplifying the relation of the Cauchy stress tensor (18), to [35, 46]

As explicated in the previous sections, the displacement field inside the in vivo tumor could be measured in a noninvasive way. With the purpose of estimating Mooney–Rivlin hyperelastic parameters of the tumor just by using the displacement quantities, we manipulate (25), that is, the stress-stretch relation of the hyperelastic model,

The uniaxial load applied to the medium in alignment with the first coordinate axis practically generates the axial strain in the same direction; therefore, the strain and relation between the stress and strain could be expressed as [47, 48]

With respect to the stress-strain relation, (26), the iterative algorithm propounded for estimating Mooney–Rivlin hyperelastic parameters of a completely unknown interior tissue (i.e., tumor) could be described in the following steps:(1)Image the tissue and its adjacent mediums by the use of clinical US imaging system before/while applying a sinusoidal load with low frequency (to annul the inertia) to the exterior medium. In other words, record the relevant US RF signals or images.(2)Extract the displacement fields inside the tissue at some sequential step times (i.e., eight consecutive instants) from the recorded US RF signals or images.(3)Simulate the tissue and its neighboring mediums by making use of one of the FEM softwares corresponding to(a)the recorded predeformation US images(b)the loading specifications(c)the boundary conditions It is assumed that just the tumor and its mechanical characteristics are unidentified. The mechanical parameters of almost all healthy soft tissues have been reported in the literature. The parameters have been predominantly estimated by performing meticulous in vivo or ex vivo experiments.(4)Consider the elastic modulus, *E*, and Poisson’s ratio, *υ*, of the simulated tumor, named elastic tumor, equal 1 Pa and 0.5, respectively.(5)Compute the displacement fields inside the particular elastic tumor at the same consecutive step times by dint of the selected FEM software.(6)Calculate the real elastic modulus of the understudy tumor, *E*_{realt}, with the help of the MATLAB^{®} software (The MathWorks, Inc., Natick, Massachusetts, USA) using(a)the axial displacement quantities of several points of the tumor at some consecutive instants (based on the achieved outcomes, the axial displacement values of twelve points of the tumor at eight step times), **Y**_{realt}(b)the axial displacement values of the identical points of the elastic tumor at the same moments, **D**(c)the relation [47–49](7)Specify a set of strains, , and compute the correspondent stresses, **σ**, using (28) in accordance with (in the first iteration):(a)The strain field in the tumor could be roughly approximated from the displacement estimates.(b)The achieved results imply the selection of a set of high strains in the first iteration since the strain values are being modified in some steps of the proposed iterative algorithm; therefore, the strain values could be reduced uniformly.(c)Slight changes in the stress values, for instance by the use of the normal distribution function, might cause the stress and strain values to conform more effectively to the Mooney–Rivlin hyperelastic model assigned to the tumor,(8)Compute hyperelastic parameters of the Mooney–Rivlin model, *C*_{10} and *C*_{01}, using the stress and strain sets and the relation between stress and strain, as represented,with the help of MATLAB algorithms, for instance, the regression algorithms.(9)Assign the estimated hyperelastic parameters to the simulated tumor and compute the displacement fields inside the tumor at the determinate step times by means of the FEM software.(10)Calculate the elastic constant of the simulated tumor, *E*_{estt}, as explained previously in step 6, using the axial displacement quantities (of the appointed points at the selected step times) of the simulated tumor, **Y**_{estt}, and the elastic tumor, **D** (determined in step 6), by applying the relation(11)Appraise the estimated hyperelastic parameters for the tumor through considering the following:(a)The error of the computed axial displacement values of the selected points at the specified moments, **Y**_{estt}, by comparing them with the correspondent displacement quantities estimated from the recorded US RF signals or images, **Y**_{realt}(b)The error of the elastic constant calculated for the tumor, *E*_{estt}, by comparing it with the real elastic modulus of the understudy tumor, *E*_{realt}, estimated in step 6 The errors of the hyperelastic parameters estimated for the tumor, by comparing them with the real hyperelastic parameters of the tumor, could not be considered because it has been assumed that the tumor is entirely obscure.(12)Alter the set of strains specified in step 7 (based on the above explanation, decrease them regularly) and repeat steps 7 to 12. By reducing the strain values steadily, the error of the estimated elastic parameter for the tumor and the error of the calculated displacement quantities in the tumor are decreasing below the defined tolerance values, as illustrated,where *e*_{elastic} and *e*_{displacement} are the tolerance values and *k* represents the number of iterations of the algorithm. By decreasing the strain values chosen with regard to the mentioned conditions, the strain and stress values successively adjust more to the Mooney–Rivlin stress-strain relationship of the tumor.

#### 3. Results

##### 3.1. Soft Tissue Simulation

The breast tissue (with the dimensions of 100 × 60 × 20 mm^{3}), simulated using the FEM software ABAQUS, has been depicted in Figure 1. The breast tissue consists of three partitions, namely, fat, fibroglandular, and tumor. The Mooney–Rivlin hyperelastic model has been applied to obtain the response of simulated breast tissue to the external sinusoidal load. The Mooney–Rivlin material constants of the named breast tissues have been presented in Table 1. Since we have utilized the elastic parameter of the tumor for estimating its hyperelastic parameters, we have additionally reported the elastic parameters of the named breast tissues in Table 1. The linear and nonlinear elastic parameters have been supposed to be constant throughout each tissue partition. This set of hyperelastic parameters has been broadly utilized in the literature to simulate the breast tissue [38, 50–54].