A Pointwise Method for Identifying Biomechanical Heterogeneity of the Human Gallbladder

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(Smythe et al., , 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(Li et al., , 2013Xiong 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.
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 = Ge Ht + M, where  Table 1.
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 = λ θ ij D i1 , h ij = λ h ij 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 , κ i m (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., κ 3 1 and κ 3 2 vanish. Hence, there are a total of 13 parameters to be determined.
The in-plane Cauchy stress components at t j are: 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.
where σ φ i1 and σ θ i1 (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 , κ 1 1 , κ 1 2 , κ 1 4 , c 2 , κ 2 1 , κ 2 2 , κ 2 4 , c 3 , and κ 3 4 , but [0.01, 3] for κ 1 3 , κ 2 3 , and κ 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.
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(Li et al., , 2012(Li et al., , 2013. The stretch components of GB No. 1 over time are shown in Figure 2B. Note that the stretch-volume profiles are patient-specific. 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.

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.
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    (10) and compared with homogeneous model (Li et al., 2013) in uniform thickness h 11 = h 21 = h 31 = 2.5 mm. 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.
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, κ 2 1 and κ 2 3 , basically agree with κ 1 and κ 3 in the homogenous model, i.e., κ 2 1 ≈ (1-2) κ 1 and κ 2 3 ≈ (1-2) κ 3 . For GB 1, 19, and 21, κ 2 1 and TABLE 4 | The first principal stresses in 10 GB samples wall estimated by using the homogenous model in Li et al. (2013)   κ 2 3 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 Frontiers in Physiology | www.frontiersin.org  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.

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.
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).

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.
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.
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 κ 1 4 , c 2 , κ 2 4 , c 3 , κ 3 3 , and κ 3 4 , 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.
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). 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 9 | Relative changes in the parameters due to varied wall thickness. 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. 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 , κ 1 2 , κ 1 3 , κ 1 4 , c 2 , κ 2 3 , are κ 2 4 can be large. Hence, the segmentation error should be reduced as much as possible to 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.
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 timeconsuming (∼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.
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.
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.