Application of micro-computer tomography and inverse finite element analysis for characterizing the visco-hyperelastic response of bulk liver tissue using indentation

In-vitro mechanical indentation experimentation is performed on bulk liver tissue of lamb to characterize its nonlinear material behaviour. The material response is characterized by a visco-hyperelastic material model by the use of 2-dimensional inverse finite element (FE) analysis. The time-dependent behaviour is characterized by the viscoelastic model represented by a 4-parameter Prony series, whereas the large deformations are modelled using the hyperelastic Neo-Hookean model. The shear response described by the initial and final shear moduli and the corresponding Prony series parameters are optimized using ANSYS with the Root Mean Square (RMS) error being the objective function. Optimized material properties are validated using experimental results obtained under different loading histories. To study the efficacy of a 2D model, a three dimensional (3D) model of the specimen is developed using Micro-CT of the specimen. The initial elastic modulus of the lamb liver obtained was found to 13.5 kPa for 5% indentation depth at a loading rate of 1 mm/sec for 1-cycle. These properties are able to predict the response at 8.33% depth and a loading rate of 5 mm/sec at multiple cycles with reasonable accuracy. The visco-hyperelastic model accurately models the large displacement as well as the time-dependent behaviour of the bulk liver tissue. Mapped meshing of the 3D FE model saves computational time and captures localized displacement in an accurate manner. The 2D axisymmetric model while predicting the force response of the bulk tissue, cannot predict the localized deformations. The visco-hyperelastic model accurately models the large displacement as well as the time-dependent behaviour of the bulk liver tissue. Mapped meshing of the 3D FE model saves computational time and captures localized displacement in an accurate manner. The 2D axisymmetric model while predicting the force response of the bulk tissue, cannot predict the localized deformations.


Introduction
Many biological tissues exhibit time dependent material response and undergo large deformations which need be characterized by viscoelastic and hyperelastic material models respectively. Tensile testing requires standalone specimens and cannot predict localized response [1]. Indentation offers a non-destructive way of testing in-situ and in-vivo specimens and is particularly advantageous for testing tissues. Many studies on biological tissues are performed on animal organs due to the lack of availability of human tissues. However, certain researchers did perform experiments on human tissues. Basdogan [2] used in-vitro dynamic testing on human liver, Liu et al. [3] used static indentation in human breast tissue and Lim et al. [4] performed in-situ experiments on human liver in both static and dyamic mode. The main objective is to distinguish between health and diseased organs. Mazza et al. [5,6] compared diseased and healthy human livers using invivo aspiration experiments while Reiter et al. [7] used various experimental techniques to study hepatic fibrosis in human livers.
Compared to human tissues, the work on animal tissues is much more widespread. Samur et al. [8] have reported that the data obtained from the tests performed on live animals contain heavy noise due to breathing and other internal processes of the body. Hence, data obtained on corpses provide much stable data. Both dynamic [9,10] using mechanical stimulators and static loading were used by researchers. Indenters of spherical [11] and cylindrical [12] shape are commonly used on account of their symmetry while damage studies are usually done using mechanical graspers [13]. Among soft tissues, brain tissue characterization is relevant for studying regions of cognitive development [14], diseases such as tumors [15] and response to various types of injury [16,17]. However, the liver is highly affected by toxicity and hence has the highest number of mechanical studies. Fu et al. [18] performed tension and compression studies on porcine livers while Jordan et al. [19] used indentation. To prevent damage to tissues, tissue constructs and other artificial materials such as soft gels and rubbers which show a similar type of mechanical response as that of tissues are studied [20,21]. The advantage with these materials is that they can be tailored to any shape for testing.
When it comes to material models describing the response of soft tissues, hyperelastic and viscoelastic time dependent models are commonly used. Martins et al. [22] provided a comparative study of various hyperelastic models applied to biological tissues. The hyperelastic models generally used are Blatz model [23], Arruda Boyce model [24], logarithmic and polynomial models [25], exponential model [26], Ogden model [27] etc. The applicability of a material model depends upon the amount of deformation that the tissue undergoes during mechanical testing. Liu et al. [28] developed properties for large deformation studies using a nonlinear viscoelastic model. Valtorta et al. [29] proposes that the assumption of homogeneity and isotropy holds good for bulky organs like liver and kidney, which have limited or no reinforcement by muscular fibres. The use of finite element (FE) analysis has been prevalent as of late, with inverse FE analysis being a tool for optimization. Some of the studies using this method were done by Hu et al. [30] on pig liver and Schwartz et al. [31] on deer liver. 3D modelling using imaging techniques such as CT scans [32] and 3D ultrasound imaging [19] are generally used for specimen that are of non-standard shape.
In this work, in-vitro indentation is performed using a spherical tip indenter on lamb liver specimens under static loading using a Dynamic Mechanical Analyzer (DMA). Inverse FE analysis is performed using a visco-hyperelastic 2D model in ANSYS software. Optimization of material parameters is done by minimizing the RMS error of the material response force between experimental and FE results. 3D models of the tissue samples are developed using an imaging technique on 2D scans obtained by Micro-CT imaging of the samples. FE analysis is performed on the 3D models using the optimized material properties to compare the shear response of the 2D and 3D models.
In this article, the sample preparation, experimentation and development of 3D models from the Micro-CT scans is explained in detail in the Sect. 2: 'Method' and the material models describing the material response, the development of 2D and 3D FE models and the optimization technique used are described under Sect. 3: 'Theory' . The 'Results' are explained in Sect. 4, followed by 'Discussion' and 'Conclusion' in Sects. 5 and 6 respectively.

Method
The in-vitro indentation experimental procedure begins with the preparation of the specimen for the test. After the specimen is prepared, displacement control tests, where in a displacement is applied and the resultant force is measured are performed.

Specimen preparation
Whole liver of a Katahdin lamb is procured within an hour of slaughter and removal from the local grocery and stored in salt solution to maintain its moisture and prevent it from drying out. Then a circular piece of specimen is cut from the whole liver using a circular cutter whose diameter is almost the same as that of the glass beaker used to hold the specimen (80 mm). The cut specimen exactly fits the beaker as shown in Fig. 1. This sample is stored in salt water solution before and in between tests.

Experimental setup for displacement control testing using DMA
The setup of the liver indentation on the DMA is shown in Fig. 2. The liver is placed in a transparent glass beaker with markings along its circumference to identify the location tested. The indenter is programmed to indent the material at a specified indentation rate up to a certain depth and hold its position at that depth for a certain amount of time. When the displacement is held constant, stress relaxation occurs in the material. The corresponding reaction force exerted by the specimen is measured by the DMA. A number of indentation rates and depths are tried out on various specimens.

Development of 3-dimensional model
The specimen tested along with the beaker and indenter is placed in an acrylic chamber and placed in a Micro-CT chamber for scanning. The specimen is rotated 360° about its vertical axis; as such, the X-rays are projected from every direction through its thickness and a scan of a layer of the specimen (slice) is obtained. The slice thickness was set at 0.24 mm. The slices thus obtained using Micro-CT are stacked up over one another using Medical Imaging software developed at CVIP Laboratory, University of Louisville to obtain the actual shape of the specimen. A detailed explanation for the development of a 3D model is beyond the scope of this paper. The completed specimen is shown in Fig. 3.

Viscoelastic and Hyperelastic material models
A linear viscoelastic constitutive material model can be specified in ANSYS using two time-dependent material functions-the shear modulus G(t) and the bulk modulus K(t). It is typically observed for most time dependent materials that the bulk modulus is not very time-dependent [33]. As such, this study assumes that the bulk modulus is constant (i.e. K(t) = K 0 ). The elastic shear modulus at time t = 0 (G 0 ) is then used to evaluate K 0 via the elastic relationship Poisson's ratio is equal 0.5 for incompressible materials. For this study, the Poisson's ratio is assumed to be 0.485. This is a pretty decent assumption given that most tissues are closer to the incompressible region. It has to be pointed out here that a value of 0.5 would lead to convergence issues of FE model.
The shear modulus is specified using a Prony series given by: where G ∞ is the shear modulus at t = ∞, N is the number of Prony series terms and G i and τ i are the Prony series coefficients and relaxation times, respectively. The value of 'N' is set as 4 for all the analyses. The relaxation times τ i are evenly spaced in log time according to the equation: where τ 0 is the base relaxation time (equivalent to that of the first Prony series term τ 1 ) and 'a' is a relaxation time multiplier. For any given model solution, these values are provided by the user; a discussion of the effect of these parameters is available elsewhere [34]. The initial value of G 0 is assumed from a similar work [35] and G ∞ is initially estimated using the equation: where F end is the maximum force, i.e. at the end of loading and F end is the force at the end of relaxation (i.e. at the end of the experiment). In ANSYS, the Prony series coefficients are provided in terms of normalized values i as: where 1 , 2 and 3 and 4 are Prony series coefficients which follow the relationship, i values sum to 1. The initial elastic modulus can be determined from the shear modulus and Poisson's ratio as: Neo-Hookean hyperelastic model is chosen to model large displacements which describes the strain energy density function as: where µ is initial shear modulus (same as G 0 ), d is material incompressibility parameter defined as d = 2/K 0 , J is the ratio of the deformed elastic volume to the reference volume of the materials defined as J = det (F ij ), I 1 is the modified invariant of the first invariant I 1 of the deformation tensor C ij defined as: I 1 = J −2∕3 I 1 .

2-dimensional axisymmetric modelling and optimization
In 2D modelling, the specimen is assumed to be flat. The 2D model used for the FE analysis is shown in Fig. 4. The specimen is assumed to be perfectly flat at the top and bottom surfaces and axisymmetric about the vertical axis. The model is built with a thickness of 60 mm to approximate the actual liver specimen (maximum thickness of 80 mm and a minimum thickness of 40 mm). The width of the model is the radius of the beaker used to hold the specimen (42.25 mm). The rigid indenter is modelled as a rigid curve (circular arc with radius of the indenter as 4.76 mm (3/16th inch)) that is translated by displacements specified on a pilot node which is constrained to move vertically only. A progressive mesh with finer elements under the indenter region is used. Contact elements are placed between the indenter and the top surface. The bottom surface of the specimen is constrained in the vertical direction (UY = 0) and the right edge is constrained in the horizontal direction. PLANE183 elements are used with Viscohyperelastic properties. The total number of substeps used for the solution is 90, 30 for the loading region and 60 for the unloading region. The initial shear modulus (G 0 = 4.040 kPa) is calculated from the initial elastic modulus (12 kPa) calculated in a previous study using the DMA [36] and the final shear modulus (G ∞ ) is calculated using Eq. (4). 1 , 2 and 3 are initially assumed to be 0.50, 0.25 and 0.125, respectively, which sets 4 = 0.125 and the parameters τ 0 , a were set at 0.3, 3.2183 respectively [37]. Least squares optimization procedure is adapted here wherein the root mean square (RMS) error between the experimental force response and finite element reaction force response is minimized. The technique was developed by the authors during a previous study and is mentioned in detail in a previously published work [34].

3-dimensional model in ANSYS
The 3D model developed from the Micro-CT images is then imported into ANSYS software and the irregularities further smoothened out. The loading and boundary conditions remain similar to the axisymmetric model -bottom layer is constrained in the vertical direction and the displacement is applied through a rigid indenter controlled by a pilot node. Additionally, the displacement around the edges is constrained to mimic the effect of the glass beaker. The 3D model obtained is mapped meshed to ensure faster convergence of the solution. The smoothened and meshed 3D model of a liver specimen is shown in Fig. 5. The indenter is made to be rigid and SOLID95 elements are used for the body of the tissue.

Results
The test chosen for initial analysis is shown in Fig. 6. The displacement is ramped to a value of 3 mm (5% of thickness) at a rate of 1 mm/sec and held constant for a period of 10 s before removal. The response during removal of displacement is not considered. After optimization, the properties obtained are shown in Table 1. The best set of material properties obtained is listed in Table 1. The fit of the FE results with the experimental results is shown in Fig. 7.
Given the goodness of the fit (RMS Error of 0.012 N), the optimized properties (Table 1) are then used to predict the response of the specimen tested at 2 different   (Fig. 8). One test was performed for 2 cycles at loading and unloading rates of 1 mm/sec with total test duration of 20 s and a maximum indentation depth of 5 mm (8.33% of thickness). The 2nd test had slower loading and unloading rates (0.5 mm/sec) for both cycles with partial unloading. The duration of this test is 25 s with a maximum indentation depth of 5 mm (8.33% of thickness).
The responses of both the tests are predicted using the optimized properties in Table 1 and compared with the experimental values. The predicted force response is shown in Fig. 9.
However, for the 2D model to suffice, the load-displacement relationship at any point on the tissue surface should be identical. The positional dependence is studied by running simulations at 3 different locations on the 3D model (with different local surfaces and thicknesses). The material properties used are obtained from the optimization of the 2D model (Table 1). The results are shown in Fig. 10.
In the next step, the response of the 3D model is studied with the indentation performed at the exact location as in the experiment and the shape of the indentation is visually compared. The indented model along with the scan of the original specimen is shown in Fig. 11.

Discussion
The optimization technique with the 2D axisymmetric visco-hyperelastic model consisting of the Prony series and Neo-Hookean model provided a reasonable fit of experimental data with FE data and is able to predict loading at unloading at different loading and unloading rates as well as different indentation depths. When a pure viscoelastic model was used, convergence issues had crept in during the FE analysis at higher deformations. This work is aimed at identifying whether the effort to develop a 3D model for FE analysis can be justified. Given that a 2D model is required   for a preliminary analysis, the need for a 3D model is moot given the computational effort required. This study details the complexities, assumptions and limitations involved in the development of a 3D model. Based on the results where the deformation and stresses are localized, it can be concluded that if the object of the study is material property determination, a 3D model for an isotropic homogeneous material is unnecessary. However, for positional dependency of the material response, especially for identification of abnormalities, non-homogeneous materials with anisotropic properties, a 3D model can be justified. The viscohyperelastic model is able to predict the peak loads with consistency which depends upon the peak elastic modulus is (13.5 kPa), but the fit in the loading and unloading regions had considerable variation especially at different loading rates and indentation depths. This variation can be attributed to the number of Prony series coefficients or a higher order hyperelastic model. The effect of a mapped mesh in the 3D model considerably reduced the simulation time and convergence issues which where hindering the analysis with a free mesh. The effect of location of the test did have a minor variation, which could become important in case of critical analyses. Thus, it can be stated that the effort required in the development of a functional 3D model of the tissue is worth the effort. When the analyses need to be extended to include anisotropy, the development of a 3D model is of utmost importance.

Conclusion
The major purpose of this study is to predict the material behaviour of bulk soft tissues of large deformation using a generalized visco-hyperelastic material model. By using an inverse FE model with an optimization routine, the material parameters of a 4-term prony series and Fig. 9 Fit of multi-cycle data with optimized material properties from initial test  a Neo-Hookean model were identified. It was hypothesized that the data from a single load-hold experiment can be optimized to predict the response for complex loading histories such as load-unload, varying strain rate and multiple loading cycles. The study on a 2D axisymmetric model has proved that it can be predicted with a reasonable amount of accuracy. The study then proceeded towards establishing a relationship between the 2D and 3D FE models for soft tissue study. The 3D model was developed using Micro-CT imaging and analysed using the material properties developed from 2D analysis. It has been observed that the 2D model underpredicted the response of the material when compared with the 3D FE model with the same properties. The deformed surface in the FE model visually compared well with the actual specimen and both show highly localized displacements and stresses (in the case of FE simulation).
The development of 3D model is relevant to studying the localized response for abnormality study, but an optimization process on a 3D model can be computationally impractical. One of the major limitations of this study has been the aging of the tissue specimens during experimentation. Achieving the same amount of hydration for the specimen between and during tests has been a challenge and might have contributed to some variation in results. The same could be attributed to variation of elastic modulus values of porcine liver when compared to published works. The liver surface is highly wavy and estimating an average thickness value for the 2D axisymmetric model could be a point of further study. A more accurate study would be to identify the abnormalities in the 3D images and model them as separate entities for FE studies. It would be interesting work to separate the abnormal regions and analyse them as standalone specimens and also in the context of the total specimen while understanding the contact mechanics between them.
Acknowledgements The authors would like to thank Dr. Hongjian Shi & Prof. Aly A. Farag of the CVIP Laboratory at University of Louisville for their help in 3D imaging and Prof. Michael J. Voor of University of Louisville for providing the necessary imaging equipment.
Funding The work was supported by University of Louisville, Louisville, USA.

Availability of data and material
The experimental data as text is made available as supplementary material.
Code availability No code is custom built for this work.

Conflicts of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.