## Abstract

Continuing advances in mechanobiology reveal more and more that many cell types, especially those responsible for establishing, maintaining, remodelling or repairing extracellular matrix, are extremely sensitive to their local mechanical environment. Indeed, it appears that they fashion the extracellular matrix so as to promote a ‘mechanical homeostasis’. A natural corollary, therefore, is that cells will try to offset complexities in geometry and applied loads with heterogeneous material properties in order to render their local environment mechanobiologically favourable. There is a pressing need, therefore, for hybrid experimental–computational methods in biomechanics that can quantify such heterogeneities. In this paper, we present an approach that combines experimental information on full-field surface geometry and deformations with a membrane-based point-wise inverse method to infer full-field mechanical properties for soft tissues that exhibit nonlinear behaviours under finite deformations. To illustrate the potential utility of this new approach, we present the first quantification of regional mechanical properties of an excised but intact gallbladder, a thin-walled, sac-like organ that plays a fundamental role in normal digestion. The gallbladder was inflated to a maximum local stretch of 120% in eight pressure increments; at each pressure pause, the entire three-dimensional surface was optically extracted, and from which the surface strains were computed. Wall stresses in each state were predicted from the deformed geometry and the applied pressure using an inverse elastostatic method. The elastic properties of the gallbladder tissue were then characterized locally using point-wise stress–strain data. The gallbladder was found to be highly heterogeneous, with drastically different stiffness between the hepatic and the serosal sides. The identified material model was validated through forward finite-element analysis; both the configurations and the local stress–strain patterns were well reproduced.

## 1. Introduction

As noted many years ago by Fung [1], the characterization of multi-axial mechanical properties of biological tissues is fundamental to understanding the associated mechanics. Since that time, we have also learned that quantifying the mechanics is, in turn, fundamental to understanding the many biological responses that are mediated through mechanotransduction [2]. An emerging foundational hypothesis of mechanobiology suggests that cells establish, maintain, remodel or repair the extracellular matrix to achieve a preferred mechanical environment, a so-called mechanical homeostasis [3]. Towards this end, it appears that these cells can offset complexities in geometry and mechanical loading with heterogeneous material properties to achieve this homeostasis. Whereas traditional approaches in experimental mechanics and biomechanics seek to simplify the data analysis by focusing on regions of homogeneity, the mechanobiology demands new methods that enable potentially heterogeneous material properties to be assessed, often given the native geometry of the tissue or organ.

In this paper, we present a new approach for characterizing heterogeneous properties of thin-walled soft biological structures having an overall saccular shape. Our approach integrates standard and panoramic digital image correlation (DIC) techniques [4–7] with a point-wise inverse stress analysis method [8–10]; the former is used to measure full-field surface deformations and the latter to compute full-field wall stresses. Although conceptually similar to prior approaches that focus on a small number of regions of interest [11], full-field approaches are obviously preferred. The essential feature of this combined approach is that, by knowing associated strain and stress fields independently at multiple equilibria, one can determine parameters within an assumed constitutive relation at each Gauss point within the computational domain, that is, point by point. This basic procedure is summarized in figure 1. The idea of point-wise identification has been proposed before [12–15], but this work is the first, to the best of our knowledge, biomechanical study within an advanced experimental setting.

For the purposes of illustration, we present results for the quantification of regional mechanical properties of a lamb gallbladder. This small, pear-shaped organ consists primarily of three layers: an innermost layer (mucosa) composed of an epithelial layer and underlying loose connective tissue, a middle layer (muscularis) composed mainly of smooth muscle, and an outer layer composed of fibrous connective tissue and an outer serous membrane, essentially an extension of the visceral peritoneum. The primary function of the gallbladder is to concentrate and release bile, via smooth muscle contraction, into the duodenum. In the human, the gallbladder is approximately 8 cm long and 4 cm in diameter and it distends to hold 30–50 ml of bile. Although an important organ, the gallbladder has only recently attracted careful biomechanical study [16,17]. Fundamental to these prior studies, however, has been the implicit assumption of material homogeneity based on incomplete experimental data.

## 2. Methods

### (a) Experimental set-up

Panoramic digital image correlation (p-DIC) is a new technique that was developed to measure full-field surface deformations on small arteries during *in vitro* testing [4–7]. The technique is much more general, however. A detailed description of its theoretical basis can be found in [4], whereas the latest version of the experimental set-up and data-processing methodology can be found in [7]. Briefly, the panoramic capability is achieved by placing the specimen coaxially within a concave conical mirror and viewing it from above. The surface of the mirror reflects the full lateral surface in a ring-like fashion without discontinuity. In addition, the curved reflecting surface magnifies the image approximately 5× (depending on the size of the specimen), which enables higher spatial resolution (albeit spatially varying) with ordinary off-the-shelf lenses. The specimen typically needs to be patterned in a random fashion to endow the surface with unique features that allow the tracking of sets of material points during deformation according to DIC area-based matching algorithms [18]. To obtain three-dimensional information, a fixed single camera captures information at different viewing angles (*n*_{view}≥4) by rotating a tilted flat mirror that is mounted at 45^{°} and placed between the conical mirror and the camera [7].

Data collection for this project required a modification of the original p-DIC set-up because of the size and the sac-like geometry of the gallbladder. Indeed, a new data-processing methodology was implemented to capture and merge the three-dimensional information gathered not only for the lateral surface, as in prior work [7], but also for the fundus. A schematic drawing of the experimental set-up is shown in figure 2. Once the gallbladder was placed coaxially within the conical mirror, it was possible to image the fundus directly from above via standard DIC and to image the body and part of the neck region via the reflections from the conical mirror via p-DIC.

To obtain three-dimensional data, a stereo arrangement was needed for both the standard and the panoramic DIC. Multiple views were obtained by mounting a single camera on a special rig that accurately positioned the camera along a given viewing direction (e.g. along a line *c* inclined to the desired stereo-angle with respect to the axis *z* of the conical mirror; figure 2) and then allowed a rotation over 360^{°} to collect the required set of polar-symmetrical stereo views [7]. In particular, the camera was mounted on a large custom dual-axis translational stage that was fixed on a kinematic mount (figure 2). The camera could thus frame the conical mirror from above, because the kinematic mount was attached to a large-area rotation stage (with a central circular opening) that was rigidly fixed via an aluminium rail frame to the optical bench. After positioning the conical mirror beneath the camera and centring it with respect to the rotational stage via the translational and tilting stages, the camera was gradually aligned to be coaxial with the conical mirror (i.e. aligned along the *z*-axis in figure 2). The outcome of this preliminary process was highly accurate given that the dotted calibration target on the conical mirror (figure 2) was used in intermediate steps to recover the position and orientation of the camera with respect to the global reference system via a direct linear transformation method [4]. The adopted calibration procedure was proved to retrieve three-dimensional positions of the target points with an accuracy of 10^{−2} mm (see the electronic supplementary material for details) within the entire measurement volume. Once the camera was set coaxial with the conical mirror, it was then translated in one direction (e.g. *y*-axis in figure 2) a distance *t*_{y} and then tilted along the perpendicular direction (*x*-axis in figure 2) at a desired stereo-angle *θ* (in this study *θ*∼4^{°}). The translation *t*_{y} imposed before titling was determined by the condition that axis *c* of the camera should cross axis *z* of the conical mirror at a point below the vertex of the conical mirror, because this condition ensured the most accurate reconstruction of the stereo data. Finally, once the test was started, the multi-axial stage/camera assembly was rotated 360^{°} about the *z*-axis to collect the multiple stereo views required for panoramic measurement at each equilibrium state of interest. Image acquisition could be continuous (i.e. collecting many views depending on the acquisition rate and the speed of rotation) or incremental (e.g. at four positions spaced by 90^{°}, namely Nk, Ek, Sk and Wk, with k being the load step number, as done herein).

### (b) Testing protocol

A lamb gallbladder was tested within a few hours of excision. The surface did not require any special microdissection because it was naturally smooth and devoid of loose connective tissue. Figure 3*a* shows the specimen as it appeared after excision and figure 3*b* shows the tested specimen after the bile was gently rinsed out. To allow the gallbladder to be self-standing, even at low pressures when mounted vertically within the conical mirror, it was cut close to the neck region along the dotted line indicated in figure 3*b* and then secured with a suture to a custom-made round polytetrafluoroethylene adapter. Figure 3*c* shows the slightly inflated gallbladder mounted on the cannula before placement within the conical mirror. A three-way valve was used to connect the sample both to a syringe, for distending the organ, and to a pressure transducer (0.07 mm Hg resolution, 1.5 mm Hg accuracy).

To correlate sequential images via DIC, the surface needs to be covered with a greyscale random pattern. Several preliminary tests on a phantom revealed the most appropriate speckle size and contrast for accurate image matching [19]. Towards this end, a soft sponge was first used to gently apply a coarse spot grey background that was subsequently over-sprayed with either black or white Indian ink to obtain a fine speckle pattern. The patterned specimen was then mounted within the conical mirror, which was slowly filled with phosphate-buffered solution (PBS). The gallbladder was left in the PBS at room temperature for about an hour to equilibrate following preparation and mounting. It was then preconditioned prior to testing, which began as a gradual distension to a pressure of *P*_{1}=1.6 mm Hg to recover from the collapsed state. The preconditioning protocol consisted of eight cycles from approximately 1.6–30 mm Hg at a rate of approximately 0.5 mm Hg s^{−1} with a recovery period of approximately 2–3 min between subsequent cycles. Because the accuracy of DIC measurements is affected by rigid body motions, it was not possible to calculate the strain at the reference pressure (stand-up configuration) with respect to the collapsed undeformed configuration. Visual inspection revealed, however, the presence of wrinkles in the region immediately adjacent to the ligature (later not included in the analysis) that disappeared with further pressurization; this observation supports the hypothesis of negligible deformation at the reference state. Four images (N1, E1, S1, W1) were captured at this first configuration by rotating the camera from 0^{°} to 90^{°} to 180^{°} and then 270^{°}. Eight additional loaded configurations were imaged at *P*_{2}=2.5, *P*_{3}=3.5, *P*_{4}=8, *P*_{5}=13, *P*_{6}=20, *P*_{7}=30, *P*_{8}=40 and *P*_{9}=50 mm Hg, and the corresponding sets of four views were recorded at each pressure. About 3 min elapsed between each configuration. After pressure loading, the speckle pattern was gently removed with PBS, and the wall thickness was measured with a micrometer (*accuracy*±0.001 mm) at 12 points of the sample (figure 4). The thickness information was not used in subsequent calculations, but was reported here only to point out differences between the hepatic and serosal sides.

### (c) Image and data processing

Three-dimensional information along the entire surface at each pressurized configuration *k* (=1,…,9) was extracted from the four (i.e. Nk, Ek, Sk and Wk) images captured at the 0^{°}, 90^{°}, 180^{°} and 270^{°} rotations of the camera about the *z*-axis. Each view was inclined relative to the opposite one at a stereo-angle 2*θ* of about 8^{°}. Note that this p-DIC stereo-angle increases with the sample size, going from 2*θ*<2^{°} for mouse arteries having a diameter of approximately 1 mm [4–7] to 2*θ*=8–12^{°} for a sample with a diameter of 20–30 mm, as in this study. The relatively large stereo-angle herein improved the accuracy of reconstruction of the radial coordinate *r*, because out-of-plane measurement accuracy (*z*-coordinate for standard DIC, *r*-coordinate for p-DIC) is a strong function of the stereo-angle [20]. Actually, when using 2*θ*≥8^{°}, the reconstruction error on the stereo plane (see [7] for details) is almost unperceivable, thus two stereo views can be sufficient for three-dimensional measurement (i.e. a dual-fixed camera arrangement could have been used in this work). Yet, while the curved surface of the conical mirror virtually increases the real stereo-angle [4], hence allowing accurate reconstruction of the body of the gallbladder, the fundus was still reconstructed via a standard DIC system with an effective stereo-angle of 8^{°}. Because this angle is too small for accurate standard DIC measurement [20], a four-view system was finally adopted to improve the accuracy of the reconstruction of the fundus by adopting a double-parallax scheme.

A typical image (e.g. Sk) contained information on both the fundus and the body of the gallbladder. The fundus region (figure 5*a*) was processed via standard procedures of stereo DIC [18] together with the homologous counterparts in images Nk, Ek and Wk to yield three-dimensional information at the fundus (figure 5*c*) at each pressurized configuration. In particular, the images were analysed with a free-shaped grid of 1668 points over an approximately 4×10^{4} pixel^{2} area. The lateral surface was imaged in a ring-like area of approximately 8×10^{5} pixel^{2} (figure 5*b*) that was first unwarped in a 1300×1100 pixel^{2} image and then processed with the p-DIC algorithms over a 120×100 point grid with a spacing of five pixels between contiguous control points (see details in [4]). This procedure reconstructed the full lateral surface as reported in figure 5*d*. A correlation subset size of 21×21 pixel^{2} over an analysis window of 41×41 pixel^{2} was used in all DIC-based procedures. The top and lateral surface data points were automatically reconstructed in the same reference frame, because both standard and panoramic DIC used the same calibration target (dotted pattern on the conical mirror) to define the world coordinate system (WCS). The entire procedure was first applied to a spherical target to evaluate metrological performances in terms of accuracy in shape reconstruction and deformation measurement as well as the error in strain calculation. The complete benchmarking of the experimental method is reported as the electronic supplementary material.

The full experimental point data (figure 5*e*) could not be used directly to build a mesh for the subsequent point-wise analysis and identification because they were obtained from two different grids, that is, the resulting mesh would present badly shaped and highly distorted elements at the fundus–body boundary. For this reason, a data smoothing and remeshing procedure was implemented in the Rhinoceros three-dimensional CAD environment to create an optimal mesh for the further calculations. In particular, the whole set of data (*x*(*i*),*y*(*i*),*z*(*i*)) was rewritten as
*tanα*(*i*)=*y*(*i*)/*x*(*i*), *D*(*i*) is the Euclidean distance of the *i*th point of the point cloud from an arbitrarily chosen ‘vertex’ *V* close to the centre of the fundus and par(*i*) is the parameter to be smoothed. This coordinate transformation allowed the ‘developed’ surface of the sac-like geometry to be fitted easily with a non-uniform rational B spline (NURBS). In particular, the smooth geometry of the undeformed configuration (*P*_{1}=1.6 mm Hg) was obtained by plotting *X*(*i*), *Y* (*i*) and *Z*(*i*) with *par*(*i*)=*D*(*i*) and then fitting the developed three-dimensional data with a NURBS of proper *U* and *V* spans and stiffness (figure 5*f*). The original raw data points were then deleted, and a new set of developed data obtained from an optimal mesh generated in a finite-element environment was projected onto the NURBS surface (figure 5*g*,*h*). The projected mesh points were finally undeveloped into the original *x*, *y*, *z* in the WCS to give the remeshed smooth undeformed configuration that was used in the subsequent calculations. The fitting error of the NURBS-based procedure is of the order of 10^{−2} mm (see electronic supplementary material for details). The final mesh contained 7096 nodes and 7156 elements.

All the other configurations were obtained analogously by starting with raw data of each deformed configuration and by smoothing the related displacement fields *u*, *v*, *w* (i.e. imposing *par*(*i*)=*u*(*i*),… etc.) following the procedure sketched in figure 5.

### (d) Inverse stress analysis

An inverse elastostatic method [9,21] was used to compute wall stresses (actually, stress resultants or tensions) at each pressurized state. This method formulates the equilibrium problem on a known deformed configuration and solves for the unloaded configuration, thereby determining the stress in the given deformed state. Theoretical underpinnings and numerical aspects of this method are found elsewhere [8,22–24]. This method applies to both solids and structures, and formulations for membranes and shells have been developed [9,21]. For thin-walled, sac-like structures, the inverse approach has a unique attribute: it can predict wall tension without invoking actual material properties [9]. This attribute arises because such structures are (nearly) statically determinate; wall stress depends primarily on the current geometry and applied loads. The inverse approach can capitalize on this attribute because it solves the equilibrium problem directly on the deformed geometry. As demonstrated previously [10,25,26], wall tension computed via such an inverse analysis is practically insensitive to constitutive functions or specific values of the material parameters used in the analysis. This remarkable feature enables us to obtain the wall tension independent of information on strain.

Because wall tension is insensitive to the constitutive relation in this inverse analysis, we can use any reasonably assumed relation to facilitate the analysis, notwithstanding the accuracy of the ensuing ‘unloaded configurations’ or ‘inverse displacements’. In this work, we used a two-dimensional neo-Hookean model described by the energy function *μ*_{1} and *μ*_{2} are material parameters having the dimension of force per unit length; *I*_{1}=*tr***C** and *I*_{2}=*det***C** are principal invariants of the surface right Cauchy–Green tensor **C**=**F**^{T}**F**. We chose *μ*_{1}=50 N mm^{−1} and *μ*_{2}=20 N mm^{−1}. These parameters are much higher than typical values of tissue stiffness, and thus the model yielded unrealistically small displacements; however, the accuracy of displacement is inconsequential to stress prediction. The gallbladder was modelled as a shell structure for which the bending and transverse shear responses were also stipulated. As done previously [10,25], we used an approximate bending relation derived from the surface energy density. The bending moment **m** was assumed to be a linear function of curvature *ρ* as **m**=**H**** ρ**, where the bending moduli

**H**=

*h*

^{2}/12 (4∂

^{2}

*w*/∂

**C**∂

**C**|

_{C=I}), which for the neo-Hookean function was worked out to be

*μ*

_{1}. Definitions of bending moment, curvature, transverse shear force and transverse shear strain can be found in [21]. Zero displacement conditions were applied at boundary nodes. Note that this neo-Hookean model was merely a vehicle for computing wall tension; it had nothing to do with the actual tissue behaviour, which was described by an exponential-type function introduced later.

### (e) Isotropy test

If the material is isotropic, the tension **S** and the surface deformation tensor **C** obey a universal relation [27], namely
**SC**–**CS** is a 2×2 skew-symmetric tensor having effectively one non-zero element, the 1–2 or (2–1) component, which is also an invariant with respect to change of basis. The right-hand side of (2.2) is the sum of normalized values of **SC**–**CS** over all pressurized states (the superscript *k* indicates states). This indicator was computed at every Gauss point.

### (f) Constitutive model

As shown later, the isotropy indictor (II) turned out to be close to zero everywhere except near the boundary, and, for that reason, the gallbladder was deemed nearly isotropic. Given that soft tissues also typically exhibit an exponential-type stiffening, we considered an exponential-type isotropic energy function
*w* is a surface density function for which the parameter *μ* has the dimension of force per unit length; the exponent *a* is dimensionless. In (2.3), *λ*_{1} and *λ*_{2} are the in-plane principal stretches, and *λ*_{3} is the thickness stretch. The material is assumed to be incompressible, thus **C**=**I**) is given by
*a*)*μ* may be used as an indicator of tissue stiffness. For this reason, this parameter was reported, in addition to *μ* and *a*.

### (g) Parameter identification

Constitutive parameters were identified point-wise, or, in actuality, element by element. We took each element as a material sample over which the local property was assumed homogeneous. The objective function *F* to be minimized consisted of tension data from four Gauss points,
*P*_{1} and *P*_{2} are the tension functions defined in equations (2.4) and (2.5), and the barred variables stand for the inversely computed principal tensions—effectively representing experimental data. Superscripts *k* and *j* indicate the pressurized state (*k*=2,…,9) and Gauss point (*j*=1,2,3,4), respectively. States 2–9 were used, thus each objection function enlisted 8×4=32 pairs of tension–strain data. The local optimization problem for two material parameters per element was solved using a Levenberg–Marquardt algorithm in Matlab. For completeness, the boundary layer was included in regression; the associated parameters were not used in subsequent verification studies, however.

### (h) Forward verification

The identified heterogeneous parameters were mapped to the reference configuration, which allowed the deformed configurations to be predicted for each pressure load using a standard forward finite-element analysis. The gallbladder was modelled as a shell structure whereby the in-plane response was described by equation (2.3), and the flexure bending was approximated using a linear moment–curvature relation with a bending modulus

## 3. Results

### (a) Strain and surface tension

Figure 6 shows multiple views of the strain field at state 9 (*P*_{9}=50 mm Hg) as determined via the combined standard and panoramic DIC; note the marked regional differences. Recall that the first configuration (*P*_{1}=1.6 mm Hg), for which the gallbladder just recovered from its collapsed state, was taken as the reference. Regional strains in other states had a similar pattern, but smaller magnitudes. The hepatic surface experienced a larger strain than the serosal surface in this unconstrained *in vitro* test. Specifically, on the upper portion of the gallbladder, the hepatic strain was about five to six times higher than the serosal strain. Figure 7 shows the first principal tension at state 9 (*P*_{9}=50 mm Hg). Compared with the strains, the tensions were more uniform except at the pole and some scattered locations in the hepatic belly. Given that the tension was more or less uniform while the strain varied, it can be inferred that the material was heterogeneous.

Near the boundary, neither the stress nor the strain was expected to be reliable. The strain measurement was less accurate because of the proximity to the boundaries (i.e. an abrupt change in grey-level distribution between a speckled sample surface and flat coloured background results in a poor correlation). Stress solutions were less reliable because of the imposition of the displacement boundary condition along the edge. The influence of the boundary can be observed, in part, in the anomalies in figures 6 and 8. Based on visual inspection of stress and strain anomalies near the boundary, we identified a boundary layer of width 6 mm above the basal plane.

### (b) Isotropic test

Based on DIC-measured strains and (inversely) computed tensions, we calculated the isotropy indicator at every Gauss point. As revealed by figure 8, the value of II was small (approx. 0.01–0.05 outside the boundary layer), indicating that the magnitude of **SC**–**CS** was about 1–5% of **SC**. Near the boundary, there was a region where the absolute value of II reached 0.18, but this could be attributed to the boundary effect. Notwithstanding uncertainties owing to the boundary effect, the tissue appeared to behave nearly isotropically. Note that the observed coaxiality between **S** and **C** was not a consequence of using an isotropic model in the inverse analysis, because **C** was from DIC measurements.

### (c) Parameter distribution

Parameters *μ* and *a* (cf. equation (2.3)) were determined element-wise and are shown in figure 9. Values of *μ* ranged from 0.004 to 0.180 N mm^{−1}, with a mean of 0.021 N mm^{−1}, whereas *a* was between 2.67 and 51, with a mean of 10.56. These results were obtained by requiring both parameters to be positive. If this condition is released, we may obtain other element-wise best-fit values owing to the presence of local minima. Nonetheless, we found that the compound parameter (1+*μ*)*a* was more or less invariant with respect to the parametric range. From the distribution of (1+*μ*)*a* in figure 9, it is evident that the hepatic and serosal surfaces had distinct stiffness characteristics; the hepatic side was about four to five times more compliant despite being thicker.

### (d) Forward verification

The constitutive model with element-wise properties was implemented in a forward finite-element analysis to predict deformed states. Figure 10 compares predicted versus DIC-measured configurations. By visual inspection, the configurations were well recovered. The nodal displacement errors are listed in table 1. Over the eight deformed states, mean values of displacement error ranged from 0.04 to 0.21 mm, which are moderate compared with the maximum displacements, which were in the range of 0.34–3.85 mm. Errors in Green strain were computed at Gauss points via

To visually compare experimentally measured and inversely computed strains and tensions, we plotted distributions in states 3, 6 and 9 as representative examples of low-, medium- and high-pressure loads. Figures 11 and 12 present the strain component *E*_{1} (the circumferential strain) and the first principal tension. Clearly, the experimental strain and tension patterns were well recovered. In particular, spots of localized high tension in the hepatic surface were captured.

### (e) Local regression

Quality of regression varied from element to element depending on the local strain range. To demonstrate this, we present three representative locations: element A, at the middle hepatic surface where the strain is large (stretches up to 125%); element B, at the lower hepatic surface where stretch ranged up to 115%; and element C on the upper serosal surface, which experienced a stretch up to 105%. Figure 13 presents fitting results for these elements. Clearly, element A produced the best fit. The other two elements showed a progressively diminishing quality as the strain became smaller.

## 4. Discussion

This study appears to be the first wherein full-field heterogeneous material properties of a gallbladder were characterized element by element using measured local strains and computed stresses. We combined standard and panoramic DIC techniques, for the first time, to quantify the deforming surface of the body and fundus through multiple pressurized states and to calculate the associated surface strains. Wall tensions were computed using an inverse elastostatic method, which enabled stress analysis independent of information on the actual material properties. The parameters in an energy function were identified element by element over a surface mesh containing approximately 7000 elements. The derived heterogeneous model was shown to be able to accurately reproduce the deformed configurations and stress–strain patterns.

The combined standard and panoramic DIC approach allowed the entire shape of the specimen to be quantified with a high spatial resolution (element size of about 0.5×0.5 mm^{2} along circumferential and meridional directions) and a reconstruction accuracy of the order of 10^{−2} mm. Previous studies adopted ellipsoidal approximations of the shape of the gallbladder *in vivo* [16,17] based on global ultrasonography measures (length of the three axes of the ellipsoid), or they assumed smoothed geometries based on coarse sets of data (5.5 mm distance between cross sections) from *in vitro* MR scanner measurements. By contrast, the present DIC-based approach enabled measurement of the actual full-field three-dimensional deformation by tracking thousands of material points over the entire surface, which in turn enabled calculation of regional variations in strain with high accuracy.

The ability to calculate wall stress (or tension) without invoking realistic constitutive law hinges on the property of static determinacy, that is, the stress is determined by geometry and load. While the gallbladder used in this study had a geometry that normally admits static determinacy (e.g. wall thickness of 0.1–0.4 mm relative to an overall radius of approx. 25 mm), this condition still needs to be verified, particularly given the trilayered structure of the wall. In our implementation, we used two material models, a neo-Hookean model for inverse stress computation and an exponential model for eventual material description. These models had drastically different stiffness characteristics: the neo-Hookean model was three to four orders of magnitude stiffer. Nevertheless, stress results based on these two models were similar (figure 12). Although not a rigorous validation, this comparison justified, in part, the stipulated static determinacy (i.e. constitutive independence).

Although the main goal of this work was to demonstrate a method, there were several things worth noting about the findings on gallbladder properties. First, the gallbladder appeared to behave isotropically, as the stress (actually tension) and strain were essentially coaxial in all eight deformed states at all points. This observation is in contrast to some recent studies where the gallbladder was assumed anisotropic [16,17]. We acknowledge that coaxial stress and strain over a limited number of stress–strain protocols does not prove isotropy. Nonetheless, the data in this study do not otherwise suggest anisotropy. That the identified model (an isotropic stored energy function) can recover the deformed configurations well also supports the isotropy assumption, at least indirectly. Also in contrast to prior studies that assumed material homogeneity [16,17], the present findings suggested a highly heterogeneous behaviour. Both *μ* and *a* varied greatly over the surface. Even when factoring out the influence of wall thickness (which varied some from region to region; figure 4), the intrinsic property, *μ*/*h*, should also remain heterogeneous. Notably, there was a sharp difference (almost an abrupt jump) in stiffness between the hepatic and the serosal surfaces. Such a jump seems reasonable. The hepatic side is attached to the liver *in vivo* and thus receives additional support, that is, the effective distending pressure would be less. It can be stipulated that the hepatic strain *in vivo* would be smaller because of the peripheral constraints. Therefore, from a structure–function perspective, the hepatic surface does not have to be as stiff structurally as the serosal surface.

Is knowledge of heterogeneous properties truly critical for capturing the mechanical behaviour? In other words, if the gallbladder was treated as homogeneous, then how well can the ensuing model recover the displacements and stress–strain patterns? To answer this question, we identified a set of best-fit parameters assuming tissue homogeneity, and then applied the homogeneous model to the forward analysis. The parameters were *μ*=0.1845 N mm^{−1} and *a*=1.62, notably different from the arithmetic averages of element-wise values. Figure 14 shows the predicted configuration and distributions of the first principal tension and the strain component *E*_{1} at the ninth state (*P*_{9}=50 mm Hg). The computed configuration deviated notably from the DIC-measured geometry. The mean displacement error was 1.05±0.38 mm, in contrast to 0.21±0.10 mm for the heterogeneous model. Recall that the maximum magnitude of displacement for this state was 3.85 mm. Strain distributions were also smeared, showing little contrast between the hepatic side and the serosal side (recall the sharp difference shown in figure 6). The high tension regions on the hepatic surface observed in experiments (cf. figures 7 and 12) were not recovered. Clearly, the homogeneous model was ineffective in terms of recovering configurations and local stress–strain patterns.

Some aspects of implementation could be fine-tuned and possibly lead to improvement. We identified a boundary layer by visual inspection of stress and strain anomalies; the selected region (6 mm) was thus somewhat arbitrary. In previous studies [13,14], we proposed that a ‘stress sensitivity’ could be used to quantify boundary effects. Specifically, we could compute distributions of tension under altered parameters (or models) in the inverse step, introduce a sensitivity index (e.g. percentage difference in stress), and use this index to stratify the boundary region (based on the premise that points further way from the boundary would be less sensitive to parameter variation). This scheme was not adopted herein because it involves a somewhat arbitrary choice of sensitivity cut-off. Rather, we used the forward verification to retrospectively check the adequacy of our boundary layer specification. With regard to the local regression, we assigned a uniform weight to all states in the objective function. One possible consequence is that low strain states are less represented, thus resulting in a poorer fitting accuracy at the low strain end, which may be a reason why the strain error was larger in the first few states (table 1). Interestingly, the stress results seemed to not correlate directly with the strain error. This was understandable, because stress depended mainly on the accuracy of deformed configuration, not the strain. We could also carry out parameter identification over Gauss points instead of elements. Based on the observation that stress–strain states were similar over all Gauss points in a typical element (figure 13), we do not expect that switching to a point-wise regression would result in any significant difference. Had the mesh been coarser, a point-based approach may have been preferred. With regard to stress computation, we did not test mesh dependence, partly because the mesh was sufficiently fine (element sizes comparable to the pixel window used in the DIC protocol). Such tests should be considered in future implementations. Lastly, the accuracy of pressure measurements can be improved. The pressure transducer had an accuracy of 1.5 mm Hg, which was close to the pressure range in the first two states. Uncertainties in the pressure readings may have also contributed to the elevated strain errors therein.

In conclusion, we presented a genuine point-wise approach of material characterization that was made possible by integrating advanced DIC technologies and inverse elastostatic analysis. Although the present method was developed to study the gallbladder *in vitro*, it can be applied to many other tissues and organs of similar structure and geometry. Indeed, one of the potentially most important applications of this method would be to characterize evolving properties, as, for example, of intracranial aneurysms or models thereof. It was suggested years ago that regional variations in aneurysmal properties could be mechanically and mechanbiologically favourable [28], but there has yet to be a detailed assessment [29]. As noted earlier, without detailed information on potentially evolving, regionally varying material properties of any tissue or organ, our understanding of both the mechanics and mechanobiology will remain limited.

## Acknowledgements

We thank Mr Yuanming Luo (University of Iowa) for his assistance in preparing finite-element models for this and a related project.

- Received February 21, 2014.
- Accepted April 17, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.