# Estimating the Relative Stiffness between a Hepatic Lesion and the Liver Parenchyma through Biomechanical Simulations of the Breathing Process

Authors: S. Martínez-Sanchis; M. J. Rupérez; E. Nadal; E. Pareja; S. Brugger; D. Borzacchiello; R. López-Andújar; C. Monserrat
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

In this paper, a method to in vivo estimate the relative stiffness between a hepatic lesion and the liver parenchyma is presented. This method is based on the finite element simulation of the deformation that the liver undergoes during the breathing process. Boundary conditions are obtained through a registration algorithm known as Coherent Point Drift (CPD), which compares the liver form in two phases of the breathing process. Finally, the relative stiffness of the tumour with respect to the liver parenchyma is calculated by means of a Genetic Algorithm, which does a blind search of this parameter. The relative stiffness together with the clinical information of the patient can be used to establish the type of hepatic lesion. The developed methodology was first applied to a test case, i.e., to a control case where the parameters were known, in order to verify its validity. After that, the method was applied to two real cases and low errors were obtained.


---

## Body

## 1. Introduction

There are four main types of cells in the liver: hepatocytes (or parenchymal cells), Kupffer cells, endothelial cells, and stellate cells, which can receive stimulations to produce fibrous tissue [1]. A fibrous liver is less deformable and has less tissue with good function. The blood passes with difficulty and, therefore, is more prone to diseases. [2]. There are different diseases that can attack the liver, such as cirrhosis and hepatitis, and these attacks can be related to the type of cells. For example, hepatocellular carcinoma (HCC), which is the most common type of liver cancer, is produced by hepatocytes, the main type of cells in the liver. Some diseases can result in other more serious illnesses, such as cirrhosis or hepatic tumour. These attacks are localized in a region, and a main goal in clinical application is delimiting these regions in situ for their treatment.Traditionally, biopsy was the most effective test to identify damage in the liver. However, nowadays, there are also analytic tests, known as liver function tests (LFTs or LFs), which are analyses of blood to extract additional information about the state of liver. There are different types of test, but among these, the liver transaminases tests, i.e., aspartate transaminase (AST or SGOT) and alanine transaminase (ALT or SGPT), are the most important since they are used as biomarkers of a liver injury. Other important biomarkers, such as NAFLD (nonalcoholic fatty liver disease) fibrosis score, FIB-4 index, and BARD score, or combination of these, are used in the clinical setting [3]. However, these analytic tests do not predict the degree of the disease and the localization of the damaged tissue, which is very important for a good diagnosis [2, 4]. The localization of the damaged tissue can be proved through medical image. Computer Tomography (CT) and Magnetic Resonance (MR) are the most commonly used medical images, especially when there is suspicion of a hepatic tumour [5].Characterizing the liver through the study of its mechanical properties can help to establish the degree of the illness. Particularly, in HCC, the stiffness of the tumour tissue is greater than the rigidity of the normal tissue, which influences the progression of the illness and the probability of normal cells transformation into malignant cells [6, 7]. Thus, the difference of stiffness between the tissue of a tumour and the tissue of the parenchyma can be used as a biomarker. Nowadays, ultrasound elastography (USE) is the most used methodology to measure tissue stiffness. USE is an ultrasound-based technique for measuring liver stiffness according to the velocity of the elastic wave propagation across the tissue [8]. Theoretically, elastography can assess the stiffness of tumours. Several authors have studied the utility of elastography for the characterization of liver tumours ([9–12] and others). For focal liver lesions, elastography is still an open investigation, because this method cannot distinguish between benign and malignant lesions, although there are studies with promising results. For example, in Ma et al. [11], the sensitivity and specificity for differentiation of malign from benign lesions were 85% and 84%, respectively. Also, in Guo et al. [12], the sensitivity and specificity were 83.3% and 77.9%, respectively. However, other studies, such as those presented by Heide et al. [9] and Frulio et al. [10], showed that this differentiation was not so significant. Therefore, the results are not conclusive, and more research or an alternative methodology is necessary.Currently, there are different USE techniques, which depend on the measured physical quantity [13, 14]. That is, the displacement parallel to the normal stress is measured in strain imaging technique. However, the shear wave speed is measured in shear wave imaging technique [15]. In general, these techniques present some limitations; for example, tissue attenuation decreases the ultrasound signal as a function of depth; for this reason, assessment of deeper tissue or organs is limited. Other limitation is that fluid and fat also attenuate the propagation of the external stimulus applied at the skin surface as what happens with FibroScan™ (Echosens, Paris, France). For this reason, USE techniques are not recommended for people with abdominal ascites (accumulation of liquid, produced generally by cirrhosis) or obesity [16, 17]. In these cases, other methods, such as image techniques (CT or MR), are required. Concretely, clinical applications of USE in liver are limited because, in general, USE is susceptible to internal source of stress, i.e., cardiac and breathing. For this reason, in the case of the liver, it is better to measure the stiffness in the right lobe (far of the heart) because the throbbing of the heart can result in erroneous measurements. Therefore, the measurement of the stiffness of tumours, which are located near the heart, is also limited [18, 19]. The last limitation is that, previously, the patients need to be coached in breathing (to stop breathing at the end of end-exhale state or end-inhale state) because the deep breathing can increase the errors in the stiffness measurements [20].On the other hand, recently, Magnetic Resonance Elastography (MRE) has appeared as a novel noninvasive imaging technique to describe the biomechanical behaviour of soft tissues in vivo, reducing the limitations of the traditional USE method. For example, MRE can be applied in patients with ascites or obesity. These patients have problems with using the traditional USE method, as indicated in [21]. The largest amount of clinical studies in MRE have been based on the evaluation of chronical illness in the parenchymal liver [22]. Thomson et al. [23] affirm that MRE has been used clinically in the liver for diagnoses and for establishing the liver stage of fibrosis, and they showed that the MRE may be able to differentiate HCC tumour grade. However, the MRE method has some limitations. The iron-overload state is the most important, because this overload may result in the MRI signal being too low for shear wave detection on a MRE sequence. Other limitation of MRE is a possible error of diagnosis of parenchymal organ fibrosis confounding the causes that produce an increase tissue stiffness [22]. For example, in the liver, the causes of elevated stiffness include acute inflammation, biliary obstruction, passive hepatic congestion from cardiac failure, and hepatic venous obstruction. Although MRE proves capable of differentiating benign from malignant tumours in a variety of organ systems, an overlap of the causes of stiffness in the liver can be a limiting factor, and prospective studies involving larger numbers of patients are required for validation.The purpose of this work is to present a novel methodology for the patient-specific estimation of the relative stiffness between a hepatic lesion and the liver parenchyma. The methodology proposed uses a finite element (FE) simulation of the liver deformation during a natural process such as human breathing. This FE model is based on two CT images of the patient’s abdomen, corresponding to two stages: end-exhale and end-inhale. The patients of this study had a hepatic tumour whose relative stiffness regarding the liver parenchyma was estimated. Genetic Algorithms were used to find this relative stiffness using a constitutive model to characterize the mechanical behaviour of these tissues. An iterative process was designed, so the simulation of the deformation of the liver in breathing was performed varying the elastic constants of the model. The shape of the tumour volume in the FE simulation was compared with the shape of the tumour volume in the experimentally deformed liver (i.e., the volume in inhalation) until the overlap of these two volumes was maximum, thus measuring the relative stiffness between tumour and liver parenchyma.The presented methodology reduces some of the limitations that the elastography techniques present. For example, in the case of USE, the attenuation of the ultrasound signal as a function of the tissue depth allows the assessment of deeper tissue or organs. In the case of the MRE, the iron-overload state, which may result in the MRI signal being too low for shear wave detection on a MRE sequence, is not a problem in this methodology.


## 2. Materials and Methods

CT images of the abdomen and pelvis of two anonymous patients, who had an internal tumour in the liver, were supplied by the Unidad de Cirugía Hepatobiliopancreática y Trasplante Hepático of Hospital Universitari i Politècnic La Fe de València (CT images used to support the findings of this study have not been made available because they are property of this hospital). Two CT images were acquired from each patient, corresponding to two different states of the liver deformation related to position of diaphragm in breathing: end-exhale and end-inhale. Each patient’s abdomen was scanned with Philips Diamond Select Brilliance CT 64-slice. The scan parameters were 120 kVp and 86 mA s. CT images of the liver were acquired in DICOM format with a size of 512x512x258 voxels, with a voxel size of 0.64x0.64x1.5 mm.The software Simpleware (version 4.2; Synopsys, Inc., Mountain View, California, USA) was used to process these two series of images in order to obtain the FE meshes of the models. First, liver and tumour were segmented and the two 3D geometrical models were generated (one for the end-inhale state and other for the end-exhale state). Afterwards, a smoothing Gaussian filter was used to obtain continuous surfaces for both models. Finally, the two FE meshes corresponding to each state were generated (Figure1). The element type used for each tissue, liver, and tumour, was the linear tetrahedron. Figure 2 shows the models and FE meshes in both states, end-exhale and end-inhale, for the same patient.Figure 1
Scheme of the generation of the model FE meshes.Figure 2
Model and mesh of liver and tumour in the end-exhale state and in the end-inhale state.
(a)
Segmentation of liver and tumour in the end-exhale state (b)
FE meshes of liver and tumour in the end-exhale state (c)
Segmentation of liver and tumour in the end-inhale state (d)
FE meshes of liver and tumour in the end-inhale stateAfter constructing the FE meshes of the models, it was necessary to define the boundary conditions of the problem and a suitable constitutive model for liver and tumour in order to perform the simulation of their deformation in breathing.


### 2.1. Boundary Conditions

FE modelling allows simulating the mechanical response of the tissues of an organ by means of the relationship between a displacement and/or a force field applied to the organ and a constitutive model of the organ tissues. For that, it is necessary to establish the boundary conditions (BC) of the problem to be tackled, in this case, the BC corresponding to the simulation of the liver deformation in breathing. To obtain the BC of the problem, we propose the use of point set registration, whose main goal is to find correspondences between two sets of points and obtain the transformation that maps one of the sets to the other. There are two types of transformation: rigid and nonrigid. The second type is used in medical registration and it allows translation, rotation, scaling (included anisotropic scaling), and skews. There are different registering algorithms, but a lot of them have a high computational complexity and are limited to rigid registrations. For this reason, in this study, the Coherent Point Drift (CPD) algorithm [24] was used. CPD considers the alignment between both sets of points as a probability density estimation problem, where the first point set (Gaussian Mixture Model centroids) is fitted to the second point set (data points). CPD forces the centroids to move coherently preserving the topological structure through a transformation function, which is defined as the initial position plus a displacement function,v:(1)TY,v=Y+vYwhere Y is the first point set (Gaussian Mixture Model centroids). In the work of Myronenko and Song [24], the mean squared distance between the corresponding points after registration was used as error measure, which rises with the deformation level. In this study, the maximum error obtained was about 2%, and the authors concluded that CPD shows robust and accurate performance, especially in cases with noise, outliers, and missing points.In this work, CPD was used to align the nodes of the external surface of the liver mesh at exhalation considered as the initial state to the external nodes of the liver surface at inhalation; thus the displacements to be applied to the nodes of the liver mesh (Dirichlet conditions) were obtained. This point set registration algorithm allowed mapping the end-exhale FE mesh to the end-inhale FE mesh, creating a displacement field that defined the transition between both states and provided the boundary conditions of the problem. The displacements obtained from the liver registration (Figure3) were used as the boundary conditions and applied to the liver at the end-exhale state for the FE simulations of the breathing process. Once the boundary conditions were established, the final location of the tumour in the FE-simulated deformed state (end-inhale) would allow the estimation of the relative stiffness between the tumour and the liver parenchyma.Figure 3
External nodes from the nondeformed state (in red) and from the deformed state (in blue) of a liver.


### 2.2. Constitutive Model

Nonlinear constitutive models, which are described by a specific strain energy function (W), are used to define the biomechanical behaviour of the liver tissue. There are a lot of types of constitutive laws to describe the mechanical behaviour of the liver in the literature, for example, polynomial hyperelastic models, exponential and logarithmic hyperelastic models, and combined models. In this study, a first-order Ogden model was chosen to represent the mechanical behaviour of the liver and the tumour based on the results from Martínez-Martínez et al. [25] and Untaroiu and Lu [26]. For this model, the strain energy potential for the Ogden model [27], WO, is defined as (2)WO=μαλ¯1α+λ¯2α+λ¯3α-3+KO2J-12where N denotes the order of the model; μi and αi stand for the material elastic parameters; λ¯1, λ¯2, and λ¯3 denote the deviatoric stretches; K0 is the initial bulk modulus; and J is the determinant of the elastic deformation gradient.One limitation of this study is that the liver is prestressed inside the body at the end-exhale state. Establishing the nondeformed state of the liver inside the body is an unsolved problem of great difficulty due to the fact that all the interactions of the liver with the rest of the organs and fluids should be known in order to apply an algorithm to determine this nondeformed state. Therefore, a first approximation, the end-exhale state, was considered as a nondeformed state in our methodology. Anisotropy is considered another important factor. Most of the biological soft tissues have properties highly anisotropic, such as the cartilage or the muscular tissue, and the liver is not an exception. However, liver tissue is considered as isotropic in a large majority of studies [28, 29], as a good approximation to model its behaviour. Therefore, the mechanical behaviour of the liver parenchyma and the tumour tissue was considered isotropic.The elastic constants of the first-order Ogden constitutive model of the tumour,μT and αT, were related to the elastic constants first-order Ogden constitutive model of the liver parenchyma, μL and αL, as (3) and (4) show.(3)μT=krelμL(4)αT=αLThe purpose of this work was to find the relative stiffness factor between tumour and liver parenchyma using the first-order Ogden model,krel. Small strains are not considered in this work. However, this model must verify that, for small strains, the material elastic parameters (μ and α) are related to the modulus of rigidity (G) as (5) shows.(5)μ1α1=2GFor this study, the elastic constants of the liver parenchyma were fixed toμ1L = 364.74 Pa and α1L = 16.19 [25]. Therefore(6)μT=krelμL=364.74krel(7)αT=αL=16.19.The bulk modulusK0 was fixed to 104 Pa based on the results measured by Hostettler et al. [30] for the bulk modulus of the human liver in vivo.


### 2.3. Iterative Process and Cost Function

In this study, a FE simulation of the deformation from end-exhale sate to the target deformation, end-inhale, state was included within an iterative optimization routine based on Genetic Algorithms, which is aimed at finding the relative stiffness that defines the biomechanical behaviour of the tumour with respect to the liver parenchyma. The computation of the optimal values of the relative stiffness is based on an iterative process that is performed in two steps: the first step consists in estimating the error committed by a selected value of the stiffness factor, and the second one consists in recomputing this value to reduce the above-mentioned error. This iterative process is carried out by minimising an error function that is dependent on this relative stiffness factor called cost function. In this work, a metaheuristic method called Genetic Algorithms (GAs) [31] was chosen for the global optimization of the relative stiffness factor of the biomechanical model based on the result from a previous research, where the performance of gradient descent algorithms was compared with the performance of several evolutionary algorithms Martínez-Martínez et al. [32]. GAs belong to the family of evolutionary computation algorithms, which are inspired in biological evolution [33, 34]. GAs mimic the natural evolution of a population by allowing solutions (parents) to reproduce new ones (children) and competing for survival in the following generations. After generations, the population progresses toward an optimal solution.


#### 2.3.1. Geometric Similarity Function

In this work, a modified version of the Geometric Similarity Function (GSF) presented by Martínez-Martínez et al. [32] was chosen as a cost function in order to find the optimal stiffness factor. This modified version is shown in (8)GSF=ln⁡1-JCMHDwhere JC stands for the Jaccard coefficient and MHD stands for the modified Hausdorff distance. JC [35] measures the overlap between two volumes A and B providing values between 0, no overlap, and 1, total overlap, through(9)JC=A∩BA∪B.On the other hand, MHD [36] measures the maximum of the minimum distances between any pair of points belonging to the border of two different overlapping volumes through(10)MHD=max⁡dB¯i,dA¯iwhere dX¯(i) is the distance between the voxel i and the closest voxel of a volume X.The idea underlying the choice of this function is the comparison of the total deformed shape of two volumes as a whole, representing the deformed shape of an organ. A previous study about several coefficients carried out by Lago et al. [37], which could be used to measure the error committed in the comparison of a deformed state of a body with a reference state of that body, proved that JC and MHD coefficients performed better than other classic comparisons. Furthermore, it was proved that a combination of both coefficients provided more information about the error committed in the comparison than using them separately. Therefore, GSF was computed as a combination of both coefficients, so when overlap is maximum (JC = 1) and Hausdorff distance is minimum (MHD = 0), GSF tends to zero. The natural logarithm of the combination of both coefficients was chosen in order to distinguish values close to zero with more precision. It is important to notice that these coefficients, Jaccard coefficient and modified Hausdorff distance, which are commonly used in medical image analysis, make comparisons of volumes. This means that the FE meshes to be compared are voxelized before the comparison. In this work, GSF provided the error in the comparison between the deformed shape of the two tumour volumes, one corresponding to the tumour in the actual deformed state (end-inhale) and the other corresponding to the tumour in simulation of this situation using the biomechanical model with the proposed elastic parameters and relative stiffness provided by the GA algorithm. The smaller the committed error, the lower the GSF value.


#### 2.3.2. Genetic Algorithm

In this work, GAs were used to find the patient-specific relative stiffness factor between tumour and liver parenchyma within certain bounds. The choice of these bounds were based on values provided in the literature. These bounds were normalized (between 0 and 1) to allow the optimization to converge faster avoiding the cost function to take skewed shapes.Firstly, GAs randomly generated an initial population in the first generation which contained 50 individual candidate solutions. GAs create a new population, with 50 individuals, from the population of the previous generation (30 individuals are created through crossover and the remaining 20 through mutation). The two of the individual solutions in the current generation that have the best GSF values are chosen as elite and replicated to the next generation. The terminating condition in this case was achieved when the cost function value of the best solution was below a given threshold.The entire method was implemented in a MATLAB script (version: MATLAB TAH 2018a, MathWorks, Inc.), which communicates with the free FE software FEBio (version 2.7).


## 2.1. Boundary Conditions

FE modelling allows simulating the mechanical response of the tissues of an organ by means of the relationship between a displacement and/or a force field applied to the organ and a constitutive model of the organ tissues. For that, it is necessary to establish the boundary conditions (BC) of the problem to be tackled, in this case, the BC corresponding to the simulation of the liver deformation in breathing. To obtain the BC of the problem, we propose the use of point set registration, whose main goal is to find correspondences between two sets of points and obtain the transformation that maps one of the sets to the other. There are two types of transformation: rigid and nonrigid. The second type is used in medical registration and it allows translation, rotation, scaling (included anisotropic scaling), and skews. There are different registering algorithms, but a lot of them have a high computational complexity and are limited to rigid registrations. For this reason, in this study, the Coherent Point Drift (CPD) algorithm [24] was used. CPD considers the alignment between both sets of points as a probability density estimation problem, where the first point set (Gaussian Mixture Model centroids) is fitted to the second point set (data points). CPD forces the centroids to move coherently preserving the topological structure through a transformation function, which is defined as the initial position plus a displacement function,v:(1)TY,v=Y+vYwhere Y is the first point set (Gaussian Mixture Model centroids). In the work of Myronenko and Song [24], the mean squared distance between the corresponding points after registration was used as error measure, which rises with the deformation level. In this study, the maximum error obtained was about 2%, and the authors concluded that CPD shows robust and accurate performance, especially in cases with noise, outliers, and missing points.In this work, CPD was used to align the nodes of the external surface of the liver mesh at exhalation considered as the initial state to the external nodes of the liver surface at inhalation; thus the displacements to be applied to the nodes of the liver mesh (Dirichlet conditions) were obtained. This point set registration algorithm allowed mapping the end-exhale FE mesh to the end-inhale FE mesh, creating a displacement field that defined the transition between both states and provided the boundary conditions of the problem. The displacements obtained from the liver registration (Figure3) were used as the boundary conditions and applied to the liver at the end-exhale state for the FE simulations of the breathing process. Once the boundary conditions were established, the final location of the tumour in the FE-simulated deformed state (end-inhale) would allow the estimation of the relative stiffness between the tumour and the liver parenchyma.Figure 3
External nodes from the nondeformed state (in red) and from the deformed state (in blue) of a liver.


## 2.2. Constitutive Model

Nonlinear constitutive models, which are described by a specific strain energy function (W), are used to define the biomechanical behaviour of the liver tissue. There are a lot of types of constitutive laws to describe the mechanical behaviour of the liver in the literature, for example, polynomial hyperelastic models, exponential and logarithmic hyperelastic models, and combined models. In this study, a first-order Ogden model was chosen to represent the mechanical behaviour of the liver and the tumour based on the results from Martínez-Martínez et al. [25] and Untaroiu and Lu [26]. For this model, the strain energy potential for the Ogden model [27], WO, is defined as (2)WO=μαλ¯1α+λ¯2α+λ¯3α-3+KO2J-12where N denotes the order of the model; μi and αi stand for the material elastic parameters; λ¯1, λ¯2, and λ¯3 denote the deviatoric stretches; K0 is the initial bulk modulus; and J is the determinant of the elastic deformation gradient.One limitation of this study is that the liver is prestressed inside the body at the end-exhale state. Establishing the nondeformed state of the liver inside the body is an unsolved problem of great difficulty due to the fact that all the interactions of the liver with the rest of the organs and fluids should be known in order to apply an algorithm to determine this nondeformed state. Therefore, a first approximation, the end-exhale state, was considered as a nondeformed state in our methodology. Anisotropy is considered another important factor. Most of the biological soft tissues have properties highly anisotropic, such as the cartilage or the muscular tissue, and the liver is not an exception. However, liver tissue is considered as isotropic in a large majority of studies [28, 29], as a good approximation to model its behaviour. Therefore, the mechanical behaviour of the liver parenchyma and the tumour tissue was considered isotropic.The elastic constants of the first-order Ogden constitutive model of the tumour,μT and αT, were related to the elastic constants first-order Ogden constitutive model of the liver parenchyma, μL and αL, as (3) and (4) show.(3)μT=krelμL(4)αT=αLThe purpose of this work was to find the relative stiffness factor between tumour and liver parenchyma using the first-order Ogden model,krel. Small strains are not considered in this work. However, this model must verify that, for small strains, the material elastic parameters (μ and α) are related to the modulus of rigidity (G) as (5) shows.(5)μ1α1=2GFor this study, the elastic constants of the liver parenchyma were fixed toμ1L = 364.74 Pa and α1L = 16.19 [25]. Therefore(6)μT=krelμL=364.74krel(7)αT=αL=16.19.The bulk modulusK0 was fixed to 104 Pa based on the results measured by Hostettler et al. [30] for the bulk modulus of the human liver in vivo.


## 2.3. Iterative Process and Cost Function

In this study, a FE simulation of the deformation from end-exhale sate to the target deformation, end-inhale, state was included within an iterative optimization routine based on Genetic Algorithms, which is aimed at finding the relative stiffness that defines the biomechanical behaviour of the tumour with respect to the liver parenchyma. The computation of the optimal values of the relative stiffness is based on an iterative process that is performed in two steps: the first step consists in estimating the error committed by a selected value of the stiffness factor, and the second one consists in recomputing this value to reduce the above-mentioned error. This iterative process is carried out by minimising an error function that is dependent on this relative stiffness factor called cost function. In this work, a metaheuristic method called Genetic Algorithms (GAs) [31] was chosen for the global optimization of the relative stiffness factor of the biomechanical model based on the result from a previous research, where the performance of gradient descent algorithms was compared with the performance of several evolutionary algorithms Martínez-Martínez et al. [32]. GAs belong to the family of evolutionary computation algorithms, which are inspired in biological evolution [33, 34]. GAs mimic the natural evolution of a population by allowing solutions (parents) to reproduce new ones (children) and competing for survival in the following generations. After generations, the population progresses toward an optimal solution.


### 2.3.1. Geometric Similarity Function

In this work, a modified version of the Geometric Similarity Function (GSF) presented by Martínez-Martínez et al. [32] was chosen as a cost function in order to find the optimal stiffness factor. This modified version is shown in (8)GSF=ln⁡1-JCMHDwhere JC stands for the Jaccard coefficient and MHD stands for the modified Hausdorff distance. JC [35] measures the overlap between two volumes A and B providing values between 0, no overlap, and 1, total overlap, through(9)JC=A∩BA∪B.On the other hand, MHD [36] measures the maximum of the minimum distances between any pair of points belonging to the border of two different overlapping volumes through(10)MHD=max⁡dB¯i,dA¯iwhere dX¯(i) is the distance between the voxel i and the closest voxel of a volume X.The idea underlying the choice of this function is the comparison of the total deformed shape of two volumes as a whole, representing the deformed shape of an organ. A previous study about several coefficients carried out by Lago et al. [37], which could be used to measure the error committed in the comparison of a deformed state of a body with a reference state of that body, proved that JC and MHD coefficients performed better than other classic comparisons. Furthermore, it was proved that a combination of both coefficients provided more information about the error committed in the comparison than using them separately. Therefore, GSF was computed as a combination of both coefficients, so when overlap is maximum (JC = 1) and Hausdorff distance is minimum (MHD = 0), GSF tends to zero. The natural logarithm of the combination of both coefficients was chosen in order to distinguish values close to zero with more precision. It is important to notice that these coefficients, Jaccard coefficient and modified Hausdorff distance, which are commonly used in medical image analysis, make comparisons of volumes. This means that the FE meshes to be compared are voxelized before the comparison. In this work, GSF provided the error in the comparison between the deformed shape of the two tumour volumes, one corresponding to the tumour in the actual deformed state (end-inhale) and the other corresponding to the tumour in simulation of this situation using the biomechanical model with the proposed elastic parameters and relative stiffness provided by the GA algorithm. The smaller the committed error, the lower the GSF value.


### 2.3.2. Genetic Algorithm

In this work, GAs were used to find the patient-specific relative stiffness factor between tumour and liver parenchyma within certain bounds. The choice of these bounds were based on values provided in the literature. These bounds were normalized (between 0 and 1) to allow the optimization to converge faster avoiding the cost function to take skewed shapes.Firstly, GAs randomly generated an initial population in the first generation which contained 50 individual candidate solutions. GAs create a new population, with 50 individuals, from the population of the previous generation (30 individuals are created through crossover and the remaining 20 through mutation). The two of the individual solutions in the current generation that have the best GSF values are chosen as elite and replicated to the next generation. The terminating condition in this case was achieved when the cost function value of the best solution was below a given threshold.The entire method was implemented in a MATLAB script (version: MATLAB TAH 2018a, MathWorks, Inc.), which communicates with the free FE software FEBio (version 2.7).


## 2.3.1. Geometric Similarity Function

In this work, a modified version of the Geometric Similarity Function (GSF) presented by Martínez-Martínez et al. [32] was chosen as a cost function in order to find the optimal stiffness factor. This modified version is shown in (8)GSF=ln⁡1-JCMHDwhere JC stands for the Jaccard coefficient and MHD stands for the modified Hausdorff distance. JC [35] measures the overlap between two volumes A and B providing values between 0, no overlap, and 1, total overlap, through(9)JC=A∩BA∪B.On the other hand, MHD [36] measures the maximum of the minimum distances between any pair of points belonging to the border of two different overlapping volumes through(10)MHD=max⁡dB¯i,dA¯iwhere dX¯(i) is the distance between the voxel i and the closest voxel of a volume X.The idea underlying the choice of this function is the comparison of the total deformed shape of two volumes as a whole, representing the deformed shape of an organ. A previous study about several coefficients carried out by Lago et al. [37], which could be used to measure the error committed in the comparison of a deformed state of a body with a reference state of that body, proved that JC and MHD coefficients performed better than other classic comparisons. Furthermore, it was proved that a combination of both coefficients provided more information about the error committed in the comparison than using them separately. Therefore, GSF was computed as a combination of both coefficients, so when overlap is maximum (JC = 1) and Hausdorff distance is minimum (MHD = 0), GSF tends to zero. The natural logarithm of the combination of both coefficients was chosen in order to distinguish values close to zero with more precision. It is important to notice that these coefficients, Jaccard coefficient and modified Hausdorff distance, which are commonly used in medical image analysis, make comparisons of volumes. This means that the FE meshes to be compared are voxelized before the comparison. In this work, GSF provided the error in the comparison between the deformed shape of the two tumour volumes, one corresponding to the tumour in the actual deformed state (end-inhale) and the other corresponding to the tumour in simulation of this situation using the biomechanical model with the proposed elastic parameters and relative stiffness provided by the GA algorithm. The smaller the committed error, the lower the GSF value.


## 2.3.2. Genetic Algorithm

In this work, GAs were used to find the patient-specific relative stiffness factor between tumour and liver parenchyma within certain bounds. The choice of these bounds were based on values provided in the literature. These bounds were normalized (between 0 and 1) to allow the optimization to converge faster avoiding the cost function to take skewed shapes.Firstly, GAs randomly generated an initial population in the first generation which contained 50 individual candidate solutions. GAs create a new population, with 50 individuals, from the population of the previous generation (30 individuals are created through crossover and the remaining 20 through mutation). The two of the individual solutions in the current generation that have the best GSF values are chosen as elite and replicated to the next generation. The terminating condition in this case was achieved when the cost function value of the best solution was below a given threshold.The entire method was implemented in a MATLAB script (version: MATLAB TAH 2018a, MathWorks, Inc.), which communicates with the free FE software FEBio (version 2.7).


## 3. Results

### 3.1. Test Case

To begin with, the feasibility of the proposed methodology was proved using a test case. In this case, the constitutive model of liver parenchyma and tumour was obtained fixing to 10 the relative stiffness between both tissues (krel = 10) and using the elastic parameters proposed in Section 2.2. GAs were used to find this factor choosing the bounds shown in Table 1 for its searching. The algorithm performed a blind search of it. The final result was obtained choosing the best GSF value among the values obtained in the iterative process.Table 1
Lower and upper bounds of the search intervals for the optimization of the relative stiffness factor of the test sample.
Case Interval Test Case 1≤krel≤ 50Figure4 shows a schematic of the methodology applied for the test case. This approach is based on getting a reference deformed state for the liver and the tumour, which is compared with each deformed state simulated through GA and evaluated by means of a cost function (GSF). The reference deformed state was obtained from the model of the liver at the end-exhale state using a reference relative stiffness krel and the liver and tumour properties provided in Section 2.2, in addition to the model the displacements obtained with CPD algorithm as boundary conditions. The shape of the tumour in the deformed state of reference was compared with the shape of the tumour in each simulation obtained for each generation of parameters. This process was iterative, and when the minimal value of GSF was obtained, the relative stiffness was selected.Figure 4
Schematic of the process in the test case.The optimal value of the relative stiffness factor obtained by GAs for the test case is shown in Table2. The test case provided an overlap greater than 99%. The relative stiffness factor obtained was 10.003, and theoretically this value should be 10, so the relative error committed in the estimation of this factor was only 0.03%. This way, the capability of the algorithm to perform the estimation of the relative stiffness factor was proved.Table 2
Optimal values of the relative stiffness between tumour and liver estimated by GAs (krel) together with Jaccard coefficient (JC), modified Hausdorff distance (MHD), and Geometric Similarity Function (GSF) values for test case.
Case k r e l JC MHD (mm) GSF Test Case 10.003 0.99983 0.5493 -9.2906


### 3.2. Real Cases

After that, the methodology was applied to two real cases. Figure5 shows a scheme of the methodology used for the real cases. In this kind of cases, the two states of the liver, end-inhale state and end-exhale state, were compared through GA to find the relative stiffness factor between tumour and liver parenchyma. In the same way as before, the final result was obtained choosing the best GSF value among the obtained values in the iterative process. For these cases, the liver at the end-inhale state was the deformed reference state and each simulation was performed from the model of the liver at the end-exhale state to which the boundary conditions obtained with the CPD algorithm were applied. These BC were the displacements of the nodes of the external surface of the liver. It is important to highlight that, in this case, the value of krel was unknown, and, therefore, the committed error could only be evaluated studying the values of JC and MHD (overlap and distance between borders) in the comparison of the two tumour shapes.Figure 5
Schematic of the process in both real samples.For this part, the bounds shown in Table3 were chosen for this search.Table 3
Lower and upper bounds of the search intervals for the optimization of the relative stiffness factor of the actual livers for the Ogden model.
Case Interval Liver 1 1≤krel≤ 50 Liver 2 1≤krel≤ 50The relative stiffness factor provided by GAs for each of the actual samples is shown in Table4.Table 4
Optimal values of the relative stiffness between tumour and liver estimated by GAs (krel) together with Jaccard coefficient (JC), modified Hausdorff distance (MHD), and Geometric Similarity Function (GSF) values for each real sample.
Case k r e l JC MHD (mm) GSF Liver 1 31.186 0.85744 1.4342 -1.5874 Liver 2 41.502 0.92189 1.2886 -2.2996For these cases, the value of the overlap (JC coefficient) was about 85.7% for Liver 1 and a value of about 92.2% was obtained for Liver 2. The Hausdorff distances were 1.43 mm for Liver 1 and 1.29 mm for Liver 2, and the relative stiffness was 31.19 and 41.502 for Liver 1 and Liver 2, respectively.


## 3.1. Test Case

To begin with, the feasibility of the proposed methodology was proved using a test case. In this case, the constitutive model of liver parenchyma and tumour was obtained fixing to 10 the relative stiffness between both tissues (krel = 10) and using the elastic parameters proposed in Section 2.2. GAs were used to find this factor choosing the bounds shown in Table 1 for its searching. The algorithm performed a blind search of it. The final result was obtained choosing the best GSF value among the values obtained in the iterative process.Table 1
Lower and upper bounds of the search intervals for the optimization of the relative stiffness factor of the test sample.
Case Interval Test Case 1≤krel≤ 50Figure4 shows a schematic of the methodology applied for the test case. This approach is based on getting a reference deformed state for the liver and the tumour, which is compared with each deformed state simulated through GA and evaluated by means of a cost function (GSF). The reference deformed state was obtained from the model of the liver at the end-exhale state using a reference relative stiffness krel and the liver and tumour properties provided in Section 2.2, in addition to the model the displacements obtained with CPD algorithm as boundary conditions. The shape of the tumour in the deformed state of reference was compared with the shape of the tumour in each simulation obtained for each generation of parameters. This process was iterative, and when the minimal value of GSF was obtained, the relative stiffness was selected.Figure 4
Schematic of the process in the test case.The optimal value of the relative stiffness factor obtained by GAs for the test case is shown in Table2. The test case provided an overlap greater than 99%. The relative stiffness factor obtained was 10.003, and theoretically this value should be 10, so the relative error committed in the estimation of this factor was only 0.03%. This way, the capability of the algorithm to perform the estimation of the relative stiffness factor was proved.Table 2
Optimal values of the relative stiffness between tumour and liver estimated by GAs (krel) together with Jaccard coefficient (JC), modified Hausdorff distance (MHD), and Geometric Similarity Function (GSF) values for test case.
Case k r e l JC MHD (mm) GSF Test Case 10.003 0.99983 0.5493 -9.2906


## 3.2. Real Cases

After that, the methodology was applied to two real cases. Figure5 shows a scheme of the methodology used for the real cases. In this kind of cases, the two states of the liver, end-inhale state and end-exhale state, were compared through GA to find the relative stiffness factor between tumour and liver parenchyma. In the same way as before, the final result was obtained choosing the best GSF value among the obtained values in the iterative process. For these cases, the liver at the end-inhale state was the deformed reference state and each simulation was performed from the model of the liver at the end-exhale state to which the boundary conditions obtained with the CPD algorithm were applied. These BC were the displacements of the nodes of the external surface of the liver. It is important to highlight that, in this case, the value of krel was unknown, and, therefore, the committed error could only be evaluated studying the values of JC and MHD (overlap and distance between borders) in the comparison of the two tumour shapes.Figure 5
Schematic of the process in both real samples.For this part, the bounds shown in Table3 were chosen for this search.Table 3
Lower and upper bounds of the search intervals for the optimization of the relative stiffness factor of the actual livers for the Ogden model.
Case Interval Liver 1 1≤krel≤ 50 Liver 2 1≤krel≤ 50The relative stiffness factor provided by GAs for each of the actual samples is shown in Table4.Table 4
Optimal values of the relative stiffness between tumour and liver estimated by GAs (krel) together with Jaccard coefficient (JC), modified Hausdorff distance (MHD), and Geometric Similarity Function (GSF) values for each real sample.
Case k r e l JC MHD (mm) GSF Liver 1 31.186 0.85744 1.4342 -1.5874 Liver 2 41.502 0.92189 1.2886 -2.2996For these cases, the value of the overlap (JC coefficient) was about 85.7% for Liver 1 and a value of about 92.2% was obtained for Liver 2. The Hausdorff distances were 1.43 mm for Liver 1 and 1.29 mm for Liver 2, and the relative stiffness was 31.19 and 41.502 for Liver 1 and Liver 2, respectively.


## 4. Discussion

The method proposed in this document has been proved by means of a test case, where the relative stiffness was known and the algorithm made a blind research of it. In this case, the error was very low, only 0.03%. Moreover, the overlap between the tumour volume in the reference state and the tumour in the simulated state was very high, JC coefficient was 99.98%, and HMD was 0.55 mm. Taking into account the fact that the voxel size was 0.64x0.64x1.5 mm, these values showed that the methodology achieves good results and can be used to estimate the relative stiffness between the tumour and the liver parenchyma. This methodology had already been applied to establish the elastic constants of a constitutive model for the cornea [38] and for the breast [39]. However, this is the first time that this methodology is applied to estimate the relative stiffness between tumour and liver parenchyma.Regarding the real cases, the overlap between the tumour in the actual deformed state and in the simulated state is smaller than for the test case, but still acceptable. As the values of the relative stiffness for the real cases were unknown, the committed error can be analysed looking at the values of JC and MHD for the test case, which guarantees an error of 0.03%. For Liver 1, the overlap obtained in the comparison of the shape of the tumour volume was of 85% with a maximum distance between borders (MHD) of 1.4 mm, and for Liver 2 the overlap obtained was 92% and the distance between borders was 1.3 mm, indicating higher errors than for the test case, but still acceptable. The difference with the test case and the difference between both results in the real cases could be due to the material model used, but it could mainly be due to the quality of the segmentation. In other words, in the test case, the reference state of the liver was obtained from the same state used for performing the simulation (the liver at end-exhale). Hence, the segmentation error did not affect the search. However, for the real cases, two different states were used, end-inhale state as the reference state, and end-exhale state as the initial state for performing the simulation. Thus, two processes of segmentation were necessary. For this reason, the error committed in the segmentation processes could have affected the comparison between the tumour in the actual deformed state and the tumour in the simulated configuration.It is important to notice that the estimated relative stiffness is not the real relative stiffness because the end-exhale state of the liver was considered as a nondeformed state. As previously commented, this is an approximation that allows us to establish a suitable parameter to measure relative stiffness that can be used as a biomarker of a hepatic lesion, since liver is prestressed inside the body. It is also important to notice that the relative stiffness was obtained because Dirichlet BC were applied to the problem. However, the Ogden parametersμL and μT could also be obtained if the forces that the liver undergoes in the breathing process would be applied, as it was performed for the mechanical characterization of the cornea in vivo [38] or for the estimation of the elastic constants of the breast [39]. However, the in vivo measurement of these forces is still a challenge due to the difficult access to the internal organs. This makes it very difficult to establish boundary conditions and loads for the simulation process of liver deformation during the breathing process. Therefore, the presented work can be considered a small step in the in vivo estimation of the patient-specific stiffness of liver and tumour. In these cases, the results about the relative stiffness can be supported with other type of test, such as analytic tests or elastography in order to determine the state of the liver parenchyma since relative stiffness is a relative parameter. Although, the latter has limitations, it can be used as a test to verify the results.Keeping in mind that there are differences of relative stiffness between patients, it is interesting to use the relative stiffness between tumour and liver parenchyma as a biomarker, since these differences can arise due to different causes: the type of tumour (metastasis or hepatocarcinoma), the ground pathology (i.e., if liver is healthy or has cirrhosis), the patient having received chemotherapy previously, or the age of the patient (if the patient is older, he/she has more fibrosis). In our case, Liver 1 had suffered metastasis, and due to chemotherapy it was showing calcifications. However, Liver 2 had suffered a portal embolization, and this liver was showing cirrhosis and hepatocarcinoma. Therefore, the different nature of the tumours and liver parenchyma is reflected in the different relative stiffness obtained. It is important to notice that some authors have studied the stiffness as biomarker of tumours, but these studies have been on breast [40], prostate [41], or ovarian cancer cells [42]. For example, Hoyt et al. [41] showed that increasing tumour aggressiveness is associated with increased collagen in the healthy tissue, and this statement was supported by Fenner et al. [40]) and Xu et al. [42]. This result shows that stiffness may be used as biomarker to evaluate the relative metastasis in different types of cancer, because the stiffness can distinguish malignant cells from benign cells and the aggressivity of the cancer cells.Finally, it is important to highlight the methodology used to establish the BC of the problem. Establishing Dirichlet BC by CPD is a novel methodology that can be very useful for other models and applications in which the access to the organs is invasive, but displacements from two different deformed states must be obtained from medical images.


## 5. Conclusion

This paper has presented a computational method to obtain the patient-specific relative stiffness between tumour and liver parenchyma; this relative stiffness can be used together with the pathological information as a biomarker of the type of tumour. The method allows the measurements of the relative stiffness of tumours regarding liver parenchyma by means of the use of medical images. Two CT images of the liver in two different states of deformation are needed to estimate this relative stiffness factor, which is computed by an iterative algorithm based on the FE simulation of the human breathing. To establish the boundary conditions of the FE problem, a registration algorithm provides the displacements to be applied to reach the end-inhale state from the end-exhale state. This methodology tries to overlap the limitations that other methods present, for example, the limitations in ultrasound elastography. The method could also be used to model the behaviour of the tumour during image-guided interventions (as biopsy or radiation therapy) improving tumour targeting accuracy, reducing the margins of the irradiated tissue, or reducing the treatment duration, if the elastic constants of the liver parenchyma are known. Furthermore, this methodology could be applied to other organs, for example, lungs or prostate.As future research, this work will expand the number of cases to be studied. The main idea is getting a database large enough to definitively verify the feasibility of the proposed methodology. Moreover, the presented methodology will be applied to other organs, such as breast and prostate.



---

*Source: 5317324-2018-10-25.xml*