# A Pointwise Method for Identifying Biomechanical Heterogeneity of the Human Gallbladder

^{1}School of Engineering, University of Glasgow, Glasgow, UK^{2}Academic Surgical Unit, Royal Hallamshire Hospital, Sheffield, UK^{3}School of Mathematics and Statistics, University of Glasgow, Glasgow, UK

Identifying the heterogeneous biomechanical property of human gallbladder (GB) walls from non-invasive measurements can have clinical significance in patient-specific modeling and acalculous biliary pain diagnosis. In this article, a pointwise method was proposed to measure the heterogeneity of ten samples of human GB during refilling. Three different points, two on the equator of GB body 90° apart and one on the apex of GB fundus, were chosen to represent the typical regions of interest. The stretches at these points were estimated from ultrasound images of the GB during the bile emptying phase based on an analytical model. The model was validated against the experimental data of a lamb GB. The material parameters at the different points were determined inversely by making use of a structure-based anisotropic constitutive model. This anisotropic model yielded much better accuracy when compared to a number of phenomenologically-based constitutive laws, as demonstrated by its significantly reduced least-square errors in stress curve fitting. The results confirmed that the human GB wall material was heterogeneous, particularly toward the apex region. Our study also suggested that non-uniform wall thickness of the GB was important in determining the material parameters, in particular, on the parameters associated with the properties of the matrix and the longitudinal fibers—the difference could be as large as 20–30% compared to that of the uniform thickness model.

## Introduction

Human gallbladder (GB) is a small pear-shaped organ that is attached to the underside of the right lobe of the liver. Its function is to store and concentrate bile produced continuously by the liver. Induced by cholecystokinin (CCK), bile can be expelled from the GB to the gut to aid the digestion of fat. Cholecystitis, often due to blockage of the cystic duct by gallstones, and acalculous biliary pain are common GB diseases that affect both women and men (Cozzolino et al., 1963; Williamson, 1988). The symptoms in acalculous biliary pain disease vary widely from discomfort to severe pain, which usually follows food intake. However, the painful symptoms remain in nearly 50% patients following gallbladder removal (Cholecystectomy) (Smythe et al., 1998, 2004). This is in part due to the lack of understanding of the underlying mechanism for acalculous biliary pain.

Interestingly, many human tissues, such as artery, breast, liver, and pancreas, can develop local disease, examples include vulnerable plaque (Baldewsing et al., 2004a; Trivedi et al., 2007), atheroma in coronary and femoral arteries (Chandran et al., 2003; Baldewsing et al., 2004b; Hamilton et al., 2005), arterial stenosis (Franquet et al., 2011), and cerebral aneurysms (Zhao et al., 2011a,b). Biomechanical properties of the diseased soft tissue are different to those of the healthy ones and are often heterogeneous. Inverse methods have been developed to identify isotropic biomechanical properties (Chandran et al., 2003; Baldewsing et al., 2004a,b; Hamilton et al., 2005; Trivedi et al., 2007; Franquet et al., 2011) in terms of the Young's modulus. In studies by Zhao, Raghavan, and Lu, pointwise inverse approaches were used to reveal the anisotropic heterogeneous biomechanical properties of cerebral aneurysms (Zhao et al., 2011a,b), ascending thoracic aneurysms (Davis et al., 2015) and murine aortas (Bersi et al., 2016) on a membrane mechanic model.

Healthy human GB wall is commonly regarded as a homogenous anisotropic non-linear material in passive state, i.e., bile refilling phase (Li et al., 2012, 2013; Xiong et al., 2013). However, recent work based on *in vitro* test of a healthy lamb GB suggested that this might not be true (Genovese et al., 2014). In addition, human acalculous biliary disease can lead to increased material heterogeneity in the GB wall. In this paper, we addressed this issue by extending the homogenous anisotropic non-linear biomechanical model for human GB wall proposed in Li et al. (2013) to a heterogeneous anisotropic case. We used an inverse pointwise method to identify the heterogeneous anisotropic property at three different points on the GB wall. The method was based on an ellipsoid membrane model and an in-house developed program using MATLAB.

## Computational Models

### Geometrical Model and Stresses under Internal Pressure

A series of ultrasonic images of acalculous human GB had been scanned in 10 min interval for 60 min during the emptying phase at the Sheffield Hallamshire Teaching Hospital. A typical example is illustrated in Figure 1, marked by the three axes *D*_{1}, *D*_{2}, and *D*_{3} (*D*_{1} ≤ *D*_{2} ≤ *D*_{3}). From these images we generated the corresponding ellipsoid models, as shown in Figure 2A, which are used to estimate the GB volume.

**Figure 2. The imaged-based ellipsoid model for GB during the refilling phase, (A)** ellipsoid model with three control points on the surface, **(B)** the stretch components estimated from the ellipsoid model for GB sample No.1 listed in Table 1.

The passive biomechanical property of GB wall exhibits in the refilling phase only, hence we will focus this process. The refilling phase is the reverse process of the emptying phase. This means that for the same volume, the refilling and the emptying phases share the same ellipsoid (Li et al., 2013). The heterogeneous anisotropic biomechanical property of human GB wall in the refilling phase will be determined inversely at points 1, 2, and 3 in a GB wall. Point 1 is an intersected point of two ellipses, one is along the equator and the other is in the longitudinal direction in a meridian plane, point 2 is also on the equator but 90° apart from point 1, and point 3 is at the apex as shown in Figure 2A. In the spherical coordinate system (*r*, ϕ, θ), the coordinates of points 1, 2, and 3 are (*D*_{1}/2, 0, π/2), (*D*_{2}/2, π/2, π/2), and (*D*_{3}/2, 0, π), respectively.

Observing that the volume of the GB model was reduced by 50% by the end of emptying (Li et al., 2013), we chose the end of the emptying configuration of the GB as the reference configuration.

We interpolated the GB model with fifteen time moments throughout the refilling phase based on GB images using the geometrical similarity (Li et al., 2011). The GB wall circumferential and longitudinal in-plane stretches at point 1 at a time instant *t*_{j} is calculated with (Ragab and Bayoumi, 1998).

where the radial displacement *u*_{r} = 0.5 (*D*_{1j} − *D*_{11}), and *j* = 1,2,3…,*N*, here *N* = 15, *D*_{11} is the length of the principal axis *D*_{1} at time *t*_{1}, *D*_{1j} is the length of the principal axis *D*_{1} at time *t*_{j}. Since point 1 is on the axis of the ellipsoid, symmetry requires that ∂*u*_{ϕ}/∂ϕ = ∂*u*_{θ}/∂θ = 0. Hence, Equation (1) can be simplified to:

The incompressibility of the GB wall means that the stretch component of GB thickness must satisfy:

Similarly, the stretch components are for point 2:

and for point 3:

where *D*_{21} and *D*_{2j} are the length of principal axis *D*_{2} at time *t*_{1} and *t*_{j}, while *D*_{31} and *D*_{3j} are the length of principal axis *D*_{3} at time *t*_{1} and *t*_{j}.

These stretch components at *t*_{j} (*j* = 1,…,*N*) and point *i* (*i* = 1,2,3) can be presented simply:

The stretches during the refilling phase are plotted against the GB volume in Figure 2B at points 1, 2, and 3 for a typical GB sample. The GB volume changed with time exponentially based on an earlier model in Li et al. (2011): *V* = *G*e^{Ht} + *M*, where *G*, *H*, and *M* are parameters determined analytically using the measured GB volume and pressure at the moments *t*_{1} and *t*_{N}.

The expressions of in-plane stress components in the GB wall during the refilling phase were the same as these in the emptying phase (Li et al., 2011) since we assumed the GB material is an elastic membrane:

where *p*_{j} is the refilling pressure at time *t*_{j}, and *F*_{θ}, *F*_{ϕ}, and *F*_{n} are the functions describing the geometry of the GB:

where *K*_{1j} = *D*_{1j}/*D*_{3j}, *K*_{2j} = *D*_{2j}/*D*_{3j}, *h*_{ij} is the GB wall thickness at point *i*, and time *t*_{j}, and ${D}_{ij}={\lambda}_{ij}^{\theta}{D}_{i1}$, ${h}_{ij}={\lambda}_{ij}^{h}{h}_{i1}$. The internal pressure *p*_{j} is given by Li et al. (2013):

where *p*_{N} is the mean final bile pressure in a GB after the refilling phase chosen to be 1,466.5 Pa (11 mmHg), *p*_{1} is the bile pressure when the refilling starts, *p*_{1} = 466.6 Pa (3 mmHg), and *t*_{N} is the total time of the refilling phase. These values are estimated from *in vivo* measurements (Li et al., 2013).

### The Constitutive Model

To determine the heterogeneous material parameters of human GB wall, the structure-based anisotropic constitutive model used in Li et al. (2013) was extended so that the material parameters are location dependent. At each point, the GB wall is assumed to be composed of homogeneous matrix and two families of fibers along the circumferential and longitudinal directions, respectively, as shown in Figure 3. The strain energy functions are:

where the parameters *c*^{i}, ${\kappa}_{m}^{i}$ (*i* = 1,2,3, *m* = 1–4) are location dependent. The total number of material property constants in Equation (10) for points 1–3 should be 15 in general. However, at point 3, there are no circumferential fibers, so the second term in Equation (10) disappears, i.e., ${\kappa}_{1}^{3}$ and ${\kappa}_{2}^{3}$ vanish. Hence, there are a total of 13 parameters to be determined.

**Figure 3. The GB wall is composed of matrix and two families of fibers, the thirteen material parameters are location dependent, changing from points 1–3**.

The in-plane Cauchy stress components at *t*_{j} are:

where ${\sigma}_{i1}^{\varphi}$ and ${\sigma}_{i1}^{\theta}$ (*i* = 1,2,3) are interpreted as the initial stresses imbedded in the GB wall, which are estimated using Equations (7–9) with the internal pressure *p*_{1}.

### Comparison with Other Constitutive Models

Several phenomenological anisotropic strain energy functions have been proposed for soft tissues. Here we did not intend to be exhaustive but will choose three commonly used strain energy functions for comparisons. These include the Fung strain energy function (Ferruzzi et al., 2011):

The Choi-Vito strain energy function (Ferruzzi et al., 2011):

and the Zhou-Fung strain energy function (Zhou and Fung, 1997):

### Inverse Estimate of the Material Parameters

The material parameters in Equations (10) or (12) or (13) or (14) are selected to minimize the objective function:

The minimization was performed using the Trust-Region-Reflective algorithm in MATLAB (More and Sorensen, 1983) which terminated when the objective function value is less than 10^{− 5}. In addition, the following RMS error (RMSE) is calculated to assess the curve-fitting quality:

It should be pointed out that the optimization process was conducted at points 1, 2, and 3 simultaneously rather than separately at each point. To secure a global minimum, the initial guesses of the parameters were chosen randomly within a suitable range, such as [0.01, 10] for *c*^{1}, ${\kappa}_{1}^{1}$, ${\kappa}_{2}^{1}$, ${\kappa}_{4}^{1}$, *c*^{2}, ${\kappa}_{1}^{2}$, ${\kappa}_{2}^{2}$, ${\kappa}_{4}^{2}$, *c*^{3}, and ${\kappa}_{4}^{3}$, but [0.01, 3] for ${\kappa}_{3}^{1}$, ${\kappa}_{3}^{2}$, and ${\kappa}_{3}^{3}$. In those ranges, the optimized material parameters did not occur at the boundaries, and the curve fitting error was in the minimum. 80 initial guesses were generated randomly, followed by 80 optimization processes. The mean property constants and curve fitting errors were chosen to be the results. The detail of initial guesses on property constants optimization is given in Section Effects of Material Heterogeneity.

### The Variable Wall Thickness

The GB wall thickness is related to the stress magnitude determined from the experimental images, as shown in Equations (7) and (8). This means that even when the pressure is the same, stresses can be different due to a varied thickness. This leads to different material parameters in the strain energy function in Equation (10). We now address the issue of the variable wall thickness of GB. A three-dimensional *in vivo* measurement of wall thickness of the GB was not yet available (Engel et al., 1980; Sanders, 1980; Prasad et al., 2008; Mohammed et al., 2010; Ugwu and Agwu, 2010); however, varying thickness of GB wall was measured *in vitro* with a digital slide caliper (Su, 2005; Khan et al., 2012). A contour of GB wall thickness is illustrated in Figure 4 (Su, 2005). It is observed that the thickness of the GB apex in the fundus increases to around 5 mm maximum and the wall of the neck is as thin as 2 mm. The ratio of the maximum thickness over the thickness of the body is 1.2.

**Figure 4. A GB wall thickness profile, showing the thickest wall at the GB apex and thinnest wall near the neck, adopted from Su (2005) with permission**.

In Khan et al. (2012), 62 GB samples were divided into three age groups; (10–20) years, (21–40) years, and (41–70) years. The thicknesses of these GBs were measured manually at the fundus, body and neck. It was identified that the maximum thickness was found on the neck, and the thinnest wall is located at the fundus. For the (41–70) years group, which coincides with the patient's age group for GB surgery in the paper, the ratio of the thickness at the fundus over the thickness of the body is 0.9. This is contrary to the finding in Su (2005). These ratios are used to examine the effect of a varying thickness.

### GB Samples

The input data for our model were the geometrical parameters based on ultrasound images, and internal pressures at the start and end of emptying/refilling phase of ten GB samples from a previous study (Li et al., 2013). These geometrical parameters and internal pressures are shown in Table 1. Additionally, the geometrical parameters and pressure profile at 15 or more moments between the start and the end of refilling phase were interpolated according to the method in Li et al. (2011) and Equation (9). A uniform wall thickness was assumed at the start of refilling phase, i.e., *h*_{11} = *h*_{21} = *h*_{31} = 2.5 mm (Li et al., 2011, 2012, 2013). The stretch components of GB No. 1 over time are shown in Figure 2B. Note that the stretch-volume profiles are patient-specific.

**Table 1. Parameters of ten human GB samples, these parameters are for one dataset rather than an average of the whole set**.

In Table 1, the ejection fraction (EF) of a GB is defined as the ratio of the difference between the initial and emptied GB volumes at 30 min after venous injection of stimulator-CCK. Ethical approval for the use of data in Figure 4 and (Li et al., 2013) were approved by the ethical committees in the hospital where the studies were conducted, and the subjects gave informed consent to these studies.

### The Solution of Inverse Problem

From Table 1, we had the ellipsoid model geometrical parameters at the beginning and the end of the refilling phase, which were the same as the end (30 min after CCK) and beginning of the emptying phase from the routine ultrasound images taken in hospital, as shown in Figure 1. As the least squares method required more scattered points than the number of parameters to be estimated, the ellipsoid model (Li et al., 2011) was interpolated over 15 or more time points for the emptying phase. The internal bile pressure is then given by Equation (9). These data were used to obtain the initial guess for the optimization process. The stretch and stress components were then computed, and the objective function and the RMSE were evaluated and compared to a given criterion of 10^{− 6}. If the criterion was not satisfied, a new guess based on the Trust-Region-Reflective algorithm would be generated, and the procedure repeated until the convergent result is reached.

## Results

### Effects of the Initial Guesses on Material Property Constants

GB No.1 shown in Table 1 was randomly chosen to identify effects of initial guesses on the repeatability of inversely determined GB wall biomechanical property constants at points 1, 2, and 3. The initial guesses of the constants were generated randomly in the ranges for search of property constants mentioned in Section Inverse Estimate of the Material Parameters by normal distribution function in MATLAB and the numbers of initial guesses were specified 10, 20, 30, …, 130, respectively. The means of the determined property constants and RMSE as well as their standard error at 95% confidence level are illustrated in Table 2. The true value of these property constants should be equal to the mean ± its standard error at 95% confidence level.

**Table 2. Mean material property constants and standard error at 95% confidence level of GB No.1 with 2.5 mm uniform initial thickness at various initial guests**.

Note that the parameters at points 1–3 determined by the least squares method based on the Trust-Region-Reflective algorithm could not be repeated from one initial guess to another due to the complexity of the inverse problem. Considering the property constants determined from statistics point of view, however, the material parameters and their standard error at 95% confidence level inversely determined remained unchanged basically, especially when the number of initial guesses was 80 or more. This suggests that the biomechanical property constant values are repeatable in a statistical sense and are globally optimum. In the following sections, the property constants are extracted with 80 initial guesses at a computational cost of around 30 min. This time consumption is much less than 3.5–7.0 h based on the approach of ABAQUS 3D FEA plus MATLAB optimization solver in Li et al. (2013).

### Effects of Material Heterogeneity

The material parameters of heterogeneity inversely determined are listed in Table 3 and compared with those from the corresponding homogenous model in Li et al. (2013). The error ε in the homogenous model reflects the error in GB volume between image observation and homogenous model prediction.

**Table 3. Heterogeneous material parameters of ten GB samples determined inversely with model Equation (10) and compared with homogeneous model (Li et al., 2013) in uniform thickness h_{11} = h_{21} = h_{31} = 2.5 mm**.

For all the GBs, the material parameter associated with the matrix in the heterogeneous model is around 10 times that of the homogenous model. For GB 3, 4, 17, 29, 37, 39, and 43, the mean values of fibers-related material parameters at points 1 and 2, ${\kappa}_{1}^{2}$ and ${\kappa}_{3}^{2}$, basically agree with κ_{1} and κ_{3} in the homogenous model, i.e., ${\kappa}_{1}^{2}$ ≈ (1–2) κ_{1} and ${\kappa}_{3}^{2}$ ≈ (1–2) κ_{3}. For GB 1, 19, and 21, ${\kappa}_{1}^{2}$ and ${\kappa}_{3}^{2}$ are different from κ_{1} and κ_{3} in the homogenous model. On one hand, the material parameters in the heterogeneous model at point 1 and point 2 are similar, implying the heterogeneity of the GB wall along the circumference is small. This is especially true for GB 37 which has *D*_{1} ≈ *D*_{2}; the points 1 and 2 share the same parameters. On the other hand, the parameters at point 3 differ substantially from the other two, suggesting a strong heterogeneity from GB body to fundus.

The first principal stresses of all the GB samples are compared in Table 4 with those predicted by the homogenous model in Li et al. (2013). These stresses are extracted at point 1 since the length of an ellipsoid major axis is the shortest through that point, resulting in the highest stress level there based on Equations (7) and (8) in the ϕ direction. It is shown that the homogenous model underestimate the stresses in the wall of GB1, 3, 4, 29, 39, and 43, and overestimates them for the remaining GBs. As a result, the relative error in the first principal stresses varies in a range of −11.4%~ +10.8% in comparison with the stresses in the homogenous model.

**Table 4. The first principal stresses in 10 GB samples wall estimated by using the homogenous model in Li et al. (2013) and the heterogeneous model Equation (10) in the present paper**.

### Comparison with Other Constitutive Models

The inversely estimated parameters are shown in Tables 5–7 for the Fung, Choi-Vito, and Zhou-Fung strain energy functions, respectively. The parameters in the Zhou-Fung model were as many as 17 in total at the three points; thus the number of time instants was increased to 30 in the optimization procedure.

**Table 5. Material parameters of Fung strain energy function Equation (14) inversely determined with uniform thickness h_{11} = h_{21} = h_{31} = 2.5 mm**.

**Table 6. Material parameters of Choi-Vito strain energy function Equation (13) inversely determined with uniform thickness h_{11} = h_{21} = h_{31} = 2.5 mm**.

**Table 7. Material parameters of Zhou-Fung strain energy function Equation (14) inversely determined with uniform thickness h_{11} = h_{21} = h_{31} = 2.5 mm**.

Our results show that even though the model parameters using these strain energy functions can also be inversely determined, the errors in stress are quite large. For instance, the mean errors are 14.8, 14.6, and 14.0% for the Fung, Choi-Vito and Zhou-Fung strain energy functions, respectively, while the structure-based model Equation (10) yields a mean error of 5.0% only (Figure 5).

**Figure 5. A comparison of errors in the least-squares stress curve fitting between the present constitutive law and existing laws proposed by Fung, Choi-Vito, and Zhou-Fung, respectively**.

### Variation of the GB Wall Thickness

We notice from Table 3 that there are some large errors ranged between 5.2 and 7.8% for the parameters estimated for five GB samples: 1, 17, 21, 37, and 39. To identify the cause of the errors, the stress-volume curves of GB 3 and 39 at points 1, 2, and 3 are shown in Figure 6. The predicted stress agrees well with the observations at points 1 and 2, but not so well at point 3. In the following, we show that this is due to the uniform wall thickness assumption used in the model.

**Figure 6. Comparison of the modeled (lines) and estimated (symbols) circumferential and longitudinal stresses with the image-based ellipsoid membrane mechanic model at points 1, 2, 3, for GB 3 (A–C)**, and GB 39 **(D–F)**.

In Section Comparison with Other Constitutive Models, the GB wall heterogeneity mainly occurs in the apex region, resulting in poor agreement in the stress, as shown in Figure 6. Therefore, we altered the GB wall thickness at the apex to examine the effect of varying thickness. First, the apex thickness was changed to 3.0 mm, based on the ratio of 1.2 found in Su (2005), and kept at 2.5 mm at points 1 and 2. The extracted pointwise mechanical properties in Equation (10) are shown in Table 8.

**Table 8. Material parameters inversely determined with model Equation (10) and variable thicknesses: h_{11} = h_{21} = 2.5 mm, and h_{31} = 3.0 mm**.

The relative changes in these 13 parameters are tabulated in Table 9. The increased thickness at the apex by 20% has a considerable effect on the material parameters, with changes up to 30%, in particular, on ${\kappa}_{4}^{1}$, *c*^{2}, ${\kappa}_{4}^{2}$, *c*^{3}, ${\kappa}_{3}^{3}$, and ${\kappa}_{4}^{3}$, which are associated with the properties of the matrix and the longitudinal fibers. This is very different to the membrane theory, in which the Young's modulus is independent of the membrane thickness (Timoshenko and Woinowsky-Krieger, 1959).

We also found that an increased *h*_{31} could lower the error in the stress between the model production and the observation. If we increase *h*_{31} to 5.0 mm, the error reduces by 2.5%. Further increase in thickness does not decrease the error much, see Figure 7. Interestingly, 5.0 mm apex thickness seems to agree with measurement in Su (2005). Finally, if *h*_{31} is reduced to be 10% thinner than *h*_{11} and *h*_{21}, 2.25 m, according to Khan et al. (2012), then the errors in the stresses are greater, as shown in Figure 7. Thus, the observation that apex thickness was thinner than the GB body in Khan et al. (2012) did not agree with the results from the cohort of GB samples used here.

**Figure 7. Effect of GB wall thickness at the apex on the error in stress, the thickness at the apex is varied to be 2.25, 2.5, 3.0, 5.0, and 7.0 mm, respectively, while the thickness at the other two points 1 and 2 is kept to be 2.5 mm**.

Note that in Khan et al. (2012) post-mortem samples from “unclaimed bodies” were used and so would not have been fresh. When left *in situ* the bile will start to break down the gallbladder wall–a process known as autolysis. Therefore, the results for wall thickness might not be reliable. Samples in Su (2005), on the other hand were obtained fresh from the operating theater and washed immediately. So we would have more faith in the results in Su (2005).

### Comparison with Animal Test

Our ellipsoid model is different from the patient specific GB geometries. One may ask if such a simple model is of any practical use. To answer this question, we compared our model prediction with the *in vitro* measurements of a lamb GB (Genovese et al., 2014). In Genovese et al. (2014), a lamb GB was harvested and inflated at a pressure up to 50 mmHg, then a series points on the GB outside surface were tracked optically, and the strain fields were estimated from fitting curves of these points. The tension/stress fields were then calculated by using the elastic membrane model and solved numerically. In our model, we only used the diameters from Genovese et al. (2014) as the two minor axis lengths *D*_{1} and *D*_{2}, respectively. The stress field, hence tension, are obtained analytically from Equation (7). The results are shown in Figure 8, where the comparisons of the second Piola-Kirchhoff surface tensions are plotted for pressure *p* = 20 and 50 mmHg (we couldn't compare the results at *p* = 3.5 mmHg, as the tension profile in Genovese et al. (2014) seems to be unrealistically large, which is possibly due to a typo in the color map scale).

**Figure 8. Comparison of the peak tension from the ellipsoid model with the in vitro experimental tension of a lamb GB (Genovese et al., 2014) at the internal pressure of 20 and 50 mmHg, respectively, two plots in the top row are from Genovese et al. (2014), with permission**.

The overall agreement is encouraging; in both our model and experiments, the highest tension is found near the GB equator, and the minimum tension occurs at the apex. The values of the predicted surface tension are also in good agreement with the experimental data except some isolated tension spots due to the rapid change in the wall curvature of the lamb GB. The predicted surface tension magnitude near the GB body/equator is in a range of 0.023–0.027 N/mm, compared with 0.03–0.04 N/mm in the experiments at 20 mmHg pressure. Likewise, the predicted tension is in a range of 0.063–0.073 N/mm, compared to the range of 0.06–0.08 N/mm in the experiments at 50 mmHg.

## Discussion

Ultrasonography is a common method for monitoring GB volume variations in daily diagnosis (Dodds et al., 1985; Portincasa et al., 2003; Ugwu and Agwu, 2010). Although a detailed 3D model is more accurate, simplified geometry models are fast and therefore frequently used in clinical assessment. When an ellipsoidal model is used to estimate GB volume based on the images scanned during emptying phase, the error of the model in GB volume is about 0.8 ± 0.1 ml. This compares better to the error of 2.1 ± 0.2 ml if using sum-of-cylinder method (Dodds et al., 1985). To assess if the simplified model could predict the stress distribution of a realistic GB sample, we also compared our model prediction with the surface tension data for a lamb GB, and the overall agreement was surprisingly good.

We comment that the segmentation error of estimating the GB diameter from a GB image is usually around 4.31–7.21% (Bocchi and Rogai, 2011). To address the effect of this error on the GB wall material parameters, we introduce a random error (or noise) of 4.31–7.21% in the major axis lengths, i.e.:

where *rand* is the inner random function in MATLAB to generate a random number in value 0–1. Then we run the inverse heterogeneous problem code with these noisy data for a number of GB sample, say No. 1, 3, 17, and 21. The parameters estimated are compared in Table 10 against those without the noise.

**Table 10. Heterogeneous material parameters of four GB samples inversely determined with Equation (10) and uniform thickness h_{11} = h_{21} = h_{31} = 2.5 mm when the error in segmentation is considered**.

When noise is considered, the error in the curve fitting increased mostly by 3.1–6.7%, some can go as high as 12.3%, in comparison with the case without the noise. The influence of segmentation error on the parameters varies from one GB to another, however, the parameters at points 1, 2 are mostly likely affected by the segmentation error, particularly, changes in *c*^{1}, ${\kappa}_{2}^{1}$, ${\kappa}_{3}^{1}$, ${\kappa}_{4}^{1}$, *c*^{2}, ${\kappa}_{3}^{2}$, are ${\kappa}_{4}^{2}$ can be large. Hence, the segmentation error should be reduced as much as possible to improve the accuracy of the inverse estimation. In future, using an automatic segmentation method with small segmentation error as introduced in Bocchi and Rogai (2011) may be the way forward.

In our previous work in Li et al. (2013), the human GB wall was considered to be a non-linear composite material of matrix and two orthogonal families of fibers in the circumferential and longitudinal directions, respectively, and the material parameters were assumed to be constant. These parameters were determined inversely in Li et al. (2013) by using the FEM software-ABAQUS with a user subroutine and a MATLAB code. However, such an inverse approach is extremely time-consuming (~7 h) and unsuitable for clinical applications. In this work, we have developed a simpler approach using analytical or simpler forward solvers, which makes it possible for clinical assessment of in GB human wall disease in real time.

In addition, we extended the previous model from homogenous membrane model in Li et al. (2013) to heterogeneous model, which has significantly improved the fitting accuracy. The heterogeneity of the GB has been confirmed in the experimental work on lamb GB (Genovese et al., 2014). The inverse estimation of the heterogeneous property constants had an error less than 7% for the ten human GB samples, and the computational time was reduced by 20 times (~30 min). Further, allowing the wall thickness variation following experiments (Su, 2005), reduced the error to be less than 4%.

One potential clinical use of the model is to assessing the GB pain. In Table 4, we compare the first principal stress with the pain score associated with the CCK venous injection. It is clear that there is a strong correlation between the magnitude of the stress and the pain score. Although given the limited sample size, the homogeneous and heterogeneous model seem to do equally well in terms of pain prediction.

The limitations of our study should also be mentioned. In the study, the stretches at the points 1–3 were determined analytically during GB emptying phase. The analytical method was based on the GB volume change from images. To our best knowledge, no speckle tracking echocardiography on GB has been reported to validate our model, unlike extensive measured on human left ventricle (Edvardsen et al., 2002; Marwick, 2006; Crosby et al., 2009; Maffessanti et al., 2009; Marwick et al., 2009; Tanaka et al., 2010; Hoit, 2011; Kleijn et al., 2011). In addition, there were also no *in vitro* passive tensile tests on the specimens harvested from the body and fundus of human GB. In future, we may be able to utilize the measured strain/stretch to validate our analytical method for stretch extraction, this will make our regional GB biomechanical property identification more accurate.

Further, we only used one *in vitro* observation to determine the reference configuration of human GBs. In reality, the size of a reference configuration may not be exactly 50% of the size of totally refilled GB. It is possible that the GB reference configuration can be estimated using GB ejection fraction (EF) in cholecystokinin-cholescintigraphy (CCK-CS) (Ozden and DiBaise, 2003) or fatty meal Cholescintigraphy (FM-CS) (Al-Muqbel et al., 2010) exanimation for GB patients' *in vivo* clinical diagnosis.

Although we have investigated the thickness variation in our model, the values we used were applicable only for a healthy GB. When human GBs suffer from diseases such as acute cholecystitis, acalculous cholecystitis and ascites (Sanders, 1980; Runner et al., 2014), the GB thickness can increase significantly. Indeed, diseased GB body wall thickness was ranged in 3–5 mm (Sanders, 1980; Mohammed et al., 2010; Runner et al., 2014). How to estimate the GB wall thickness change in disease will be an important challenge for future studies.

## Conclusions

The heterogeneity of ten samples of human GB is investigated theoretically in refilling phase using a structure-based constitutive model, ellipsoidal GB and membrane in-plane mechanic model. Three different points, two on the equator of GB body with 90° apart and one on the apex of GB fundus, are chosen to evaluate the variation of the material properties. The stretches at these points are tracked over time from the routine ultrasonic images scanned at the Sheffield Hallamshire Hospital during the emptying phase. The material parameters at the three different points are determined inversely using a MATLAB code. The human GBs are found to exhibit heterogeneity, especially from GB body to its apex region. It is found that using a homogeneous model underestimate the peak stresses in the GB wall, and that a strong heterogeneity occurs from GB body to fundus. Finally, our model results indicate that the GB wall is much thicker at the apex, which clarify the contrary findings reported in the literature.

## Author Contributions

WL developed math model, designed numerical method, composed and validated code, carried out case studies, analyzed the results and drafted the manuscript. NB checked medical issues in human gallbladder, revised the draft. XL checked the math model, method, result explanation and conducted English text editing.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

The project was supported by EPSRC (Grant no. EP/G015651 EP/G028257, and EP/N014642/1).

## References

Al-Muqbel, K. M., Bani Hani, M. N., Elheis, M. A., and Al-Omari, M. H. (2010). Reproducibility of gallbladder ejection fraction measured by fatty meal cholescintigraphy. *Nucl. Med. Mol. Imaging* 44, 246–251. doi: 10.1007/s13139-010-0046-8

Baldewsing, R. A., de Korte, C. L., Schaar, J. A., Mastik, F., and van der Steen, A. F. (2004a). Finite element modelling and intravascular ultrasound elastrography of vulnerable plaques: parameters variation. *Ultrasonics* 42, 723–729. doi: 10.1016/j.ultras.2003.11.017

Baldewsing, R. A., de Korte, C. L., Schaar, J. A., Mastik, F., and van der Steen, A. F. (2004b). Finite element model for performing intravascular ultrasound elastography of human atherosclerotic coronary arteries. *Ultrasound Med. Biol.* 30, 803–813. doi: 10.1016/j.ultrasmedbio.2004.04.005

Bersi, M. R., Bellini, C., Di Achille, P., Gumphrey, J. D., Genovese, K., and Avril, S. (2016). Novel methodology for characterizing regional variations in the material properties of murine aortas. *J. Biomech. Eng.* 138, 071005. doi: 10.1115/1.4033674

Bocchi, L., and Rogai (2011). *Segmentation of Ultrasound Breast Images: Optimization of Algorithm Parameters*. Torino: Applications of Evolutionary Computation-LNCS 6624.

Chandran, K. B., Mun, J. H., Choi, K. K., Chen, J. S., Hamilton, A., Nagaraj, A., et al. (2003). A method for *in vivo* analysis for regional arterial wall property alternations with atherosclerosis: preliminary results. *Med. Eng. Phys.* 25, 289–298. doi: 10.1016/S1350-4533(02)00224-2

Cozzolino, H. J., Goldstein, F., and Green, R. R. (1963). The cystic duct syndrome. *JAMA* 185, 920–104. doi: 10.1001/jama.1963.03060120030017

Crosby, J., Amundsen, B. H., Hergum, T., Remme, E. W., Langeland, S., and Torp, H. (2009). 3-D speckle tracking for assessment of regional left ventricular function. *Ultrasound Med. Biol.* 36, 458–471. doi: 10.1016/j.ultrasmedbio.2008.09.011

Davis, F. M., Luo, Y., Avril, A., Duprey, A., and Lu, J. (2015). Pointwise characterization of the elastic properties of planar soft tissues; application to ascending thoracic aneurysms. *Biomech. Model. Mechanobiol.* 14, 967–978. doi: 10.1007/s10237-014-0646-9

Dodds, W. J., Groh, W. J., Darweesh, R. M. A., Lawson, T. L., Kishk, S. M. A., and Kern M. K. (1985). Sonographic measurement of gallbladder volume. *Am. J. Roentgenol.* 145, 1009–1011. doi: 10.2214/ajr.145.5.1009

Edvardsen, T., Gerber, B. L., Garot, J., Bluemke, D. A., Lima, J. A. C., and Smiseth, O. A. (2002). Quantitative assessment of intrinsic regional myocardial deformation by Doppler strain rate echocardiography in humans. *Circulation* 106, 50–56. doi: 10.1161/01.CIR.0000019907.77526.75

Engel, J. M., Deitch, E. A., and Sikkema, W. (1980). Gallbladder wall thickness: sonographic accuracy and relation to disease. *Am. J. Roentgenol.* 134, 907–909. doi: 10.2214/ajr.134.5.907

Ferruzzi, J., Vorp, D. A., and Humphrey, J. D. (2011). On constitutive descriptors of the biaxial mechanical behaviour of human abdominal aorta and aneurysms. *J. R. Soc. Interface* 8, 435–450. doi: 10.1098/rsif.2010.0299

Franquet, A., Avril, S., Le Riche, R., and Badel, P. (2011). Identification of heterogeneous elastic properties in stenosed arteries; a numerical plane strain study. *Comput. Methods Biomech. Biomed. Engin.* 15, 49–58. doi: 10.1080/10255842.2010.547192

Genovese, K., Casaletto, L., Humphrey, J. D., and Lu, J. (2014). Digital image correlation-based point-wise inverse characterization of heterogeneous material properties of gallbladder *in vitro*. *Pro. R. Soc. A.* 470:20140152. doi: 10.1098/rspa.2014.0152

Hamilton, A. J., Kim, H., Nagaraj, A., Mun, J. H., Yan, L. L., Roth, S. I., et al. (2005). Regional material property alterations in porcine femoral with atheroma development. *J. Biomech.* 38, 2354–2364. doi: 10.1016/j.jbiomech.2004.10.018

Hoit, B. D. (2011). Strain and strain rate echocardiography and coronary artery disease. *Circulation* 4, 179–190. doi: 10.1161/circimaging.110.959817

Khan, L. F., Naushaba, H., Paul, U. K., Banik, S., Ahmed, M., and Al-Safri, M. A. (2012). Gross and histomorphological study of thickness of the gallbladder wall. *J. Dhaka Natl. Med. Coll. Hosp.* 18, 34–38. doi: 10.3329/jdnmch.v18i1.12238

Kleijn, S. A., Aly, M. F. A., Terwee, C. B., van Rossum, A. C., and Kamp, O. (2011). Three-dimensional speckle tracking echocardiography for automatic assessment of global and regional left ventricular function based on area strain. *J. Am. Soc. Echocardiogr.* 24, 314–321. doi: 10.1016/j.echo.2011.01.014

Li, W. G., Hill, N. A., Ogden, R. W., Smythe, A., Majeed, A. W., Bird, N., et al. (2013). Anisotropic behaviour of human gallbladder walls. *J. Mech. Behav. Biomed. Mater.* 20, 363–375. doi: 10.1016/j.jmbbm.2013.02.015

Li, W. G., Luo, X. Y., Hill, N. A., Ogden, R. W., Smythe, A., Majeed, A. W., et al. (2011). A mechanical model for CCK-induced acalculous gallbladder pain. *Ann. Biomed. Eng.* 39, 786–800. doi: 10.1007/s10439-010-0205-1

Li, W. G., Luo, X. Y., Hill, N. A., Ogden, R. W., Smythe, A., Majeed, A. W., et al. (2012). A quasi-nonlinear analysis on anisotropic behaviour of human gallbladder wall. *J. Biomech. Eng.* 134:0101009. doi: 10.1115/1.4007633

Maffessanti, F., Nesser, H. J., Weinert, L., Steringer-Mascherbauer, R., Niel, J., Gorissen, W., et al. (2009). Quantitative evaluation of regional left ventricular function using three-dimensional speckle tracking echocardiography in patients with and without heart disease. *Am. J. Cardiol.* 104, 1755–1762. doi: 10.1016/j.amjcard.2009.07.060

Marwick, T. H. (2006). Measurement of strain and strain rate by echocardiography. *J. Am. Coll. Cardiol.* 47, 1313–1327. doi: 10.1016/j.jacc.2005.11.063

Marwick, T. H., Leano, R. L., Brown, J., Sun, J. P., Hoffmann, R., Lysyansky, P., et al. (2009). Myocardial strain measurement with 2-dimensional speckle-tracking echocardiography. *J. Am. Coll. Cardiol.* 2, 80–84. doi: 10.1016/j.jcmg.2007.12.007

Mohammed, S., Tahir, A., Ahidjo, A., Mustapha, Z., Franza, O., Okoye, I., et al. (2010). Sonographic gallbladder wall thickness in normal adult population in Nigeria. *South Afr. J. Radiol.* 14, 84–87. doi: 10.4102/sajr.v14i4.450

More, J. J., and Sorensen, D. C. (1983). Computing a trust region step. *SIAM J. Sci. Stat. Comput.* 4, 553–572. doi: 10.1137/0904038

Ozden, N., and DiBaise, J. K. (2003). Gallbladder ejection fraction and symptom outcome in patients with acalculous biliary-like pain. *Dig. Dis. Sci.* 48, 890–897. doi: 10.1023/A:1023039310574

Portincasa, P., Moschetta, A., Colecchia, A., Festi, D., and Palasciano, G. (2003). Measurements of gallbladder motor function by ultrasonography: towards standardization. *Digest. Liver Dis.* 35, S56–S61. doi: 10.1016/S1590-8658(03)00096-3

Prasad, M. N., Brown, M. S., Ni, C., Douek, M., Raman, S., Lu, D., et al. (2008). Three-dimensional mapping of gallbladder wall thickness on computed tomography using Laplace's equation. *Acad. Radiol.* 15, 1075–1081. doi: 10.1016/j.acra.2008.02.006

Ragab, A. R., and Bayoumi, S. E. (1998). *Engineering Solid Mechanics: Fundamentals and Applications*. Boca Raton, FL: CRC Press.

Runner, G. J., Corwin, M. T., Siewert, B., and Eisenberg, R. L. (2014). Gallbladder wall thickening. *Am. J. Roentgenol.* 202, W1–W12. doi: 10.2214/AJR.12.10386

Sanders, R. (1980). The significance of sonographic gallbladder wall thickening. *J. Clin. Ultrasound* 8, 143–146. doi: 10.1002/jcu.1870080210

Smythe, A., Ahmed, R., Fitzhenry, M., Johnson, A. G., and Majeed, A. W. (2004). Bethanechol provocation testing does not predict symptom relief after cholecystectomy for acalculous biliary pain. *Digest. Liver Dis.* 36, 682–686. doi: 10.1016/j.dld.2004.05.009

Smythe, A., Majeed, A., Fizhenry, M., and Johnson, A. G. (1998). A requiem for the cholecystokinin provocation test? *Gut* 43, 571–574. doi: 10.1136/gut.43.4.571

Su, Y. (2005). *The Mechanical Properties of Human Gallbladder*. Thesis of BEng in Mechanical Engineering, University of Sheffield, Sheffield, UK.

Tanaka, H., Hara, H., Saba, S., and Gorcsan, J. III. (2010). Usefulness of three-dimensional speckle tracking strain to quantify dyssynchrony and the site of latest mechanical activation. *Am. J. Cardiol.* 105, 235–242. doi: 10.1016/j.amjcard.2009.09.010

Timoshenko, S., and Woinowsky-Krieger, S. (1959). *Theory of Plates and Shells, 2nd Edn.* New York, NY: McGraw-Hill Book Company.

Trivedi, R. A., Li, Z. Y., U-King-Im, J., Graves, M. J., Kirkpatrick, P. J., and Gillard, J. H. (2007). Identifying vulnerable carotid plaques *in vivo* using high resolution magnetic resonance imaging-based finite element analysis. *J. Neurosurg.* 107, 536–542. doi: 10.3171/JNS-07/09/0536

Ugwu, A. C., and Agwu, K. K. (2010). Ultrasound quantification of gallbladder volume to establish baseline contraction indices in healthy adults: a pilot study. *South Afr. Radiogr.* 48, 9–12.

Williamson, R. (1988). Acalculous disease of the gall bladder. *Gut* 29, 860–872. doi: 10.1136/gut.29.6.860

Xiong, L., Chui, C. K., and Teo, C. L. (2013). Reality based modelling and simulation of gallbladder shape deformation using variational methods. *Int. J. Comput. Assist. Radiol. Surg.* 8, 857–865. doi: 10.1007/s11548-013-0821-y

Zhao, X., Raghavan, M. L., and Lu, J. (2011a). Identifying heterogeneous anisotropic properties in cerebral aneurysms; a pointwise approach. *Biomech. Model. Mechanobiol.* 10, 177–189. doi: 10.1007/s10237-010-0225-7

Zhao, X., Raghavan, M. L., and Lu, J. (2011b). Characterizing heterogeneous anisotropic properties of cerebral aneurysms with unknown stress-free geometry: a precursor to *in vivo* identification. *J. Biomech. Eng.* 133, 051008. doi: 10.1115/1.4003872

Keywords: gallbladder, strain energy function, heterogeneity, anisotropic property, constitutive law, optimization, inverse problem

Citation: Li W, Bird NC and Luo X (2017) A Pointwise Method for Identifying Biomechanical Heterogeneity of the Human Gallbladder. *Front. Physiol.* 8:176. doi: 10.3389/fphys.2017.00176

Received: 10 November 2016; Accepted: 07 March 2017;

Published: 31 March 2017.

Edited by:

Vanessa Diaz, University College London, UKReviewed by:

Haoxiang Luo, Vanderbilt University, USAAndrew James Narracott, University of Sheffield, UK

Copyright © 2017 Li, Bird and Luo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Wenguang Li, wenguang.li@glasgow.ac.uk