Volumetric quantitative optical coherence elastography with an iterative inversion method

It is widely accepted that accurate mechanical properties of three-dimensional soft tissues and cellular samples are not available on the microscale. Current methods based on optical coherence elastography can measure displacements at the necessary resolution, and over the volumes required for this task. However, in converting this data to maps of elastic properties, they often impose assumptions regarding homogeneity in stress or elastic properties that are violated in most realistic scenarios. Here, we introduce novel, rigorous, and computationally efficient inverse problem techniques that do not make these assumptions, to realize quantitative volumetric elasticity imaging on the microscale. Specifically, we iteratively solve the three-dimensional elasticity inverse problem using displacement maps obtained from compression optical coherence elastography. This is made computationally feasible with adaptive mesh refinement and domain decomposition methods. By employing a transparent, compliant surface layer with known shear modulus as a reference for the measurement, absolute shear modulus values are produced within a millimeter-scale sample volume. We demonstrate the method on phantoms, on a breast cancer sample ex vivo, and on human skin in vivo. Quantitative elastography on this length scale will find wide application in cell biology, tissue engineering and medicine. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

of disease [6]. Elastography is also rapidly becoming a part of magnetic resonance imaging instruments [7]. However, there is a relative dearth of methods that can image the mechanics of aggregates of cells or regions of tissue, in the range from tens of micrometers to tens of millimeters. Such requirements, which are necessary to migrate research in cell mechanics to a 3D in vivo or in situ context and reconcile it with whole tissue pathophysiology, exceed the range of most cellular-scale methods, and generally exceed the resolution capacity of nonoptical methods.
Optical elastography methods have recently made great strides in mapping the microscale mechanical properties of tissue [8]. Amongst these methods, optical coherence elastography (OCE), in particular, has experienced great advances [9] and is intrinsically suited to mapping mechanical properties over this intermediate range. Quantitative elastography methods, i.e., those seeking to retrieve an absolute mechanical property, typically operate by imparting a load onto the sample, and relating the resulting displacement to an elastic modulus by solving an inverse problem [10]. The inversion from displacement to elastic modulus should be consistent with physical principles that govern the deformation of a three-dimensional elastic continuum with spatially varying mechanical properties, but most current approaches in OCE do not comply with this requirement [9,11]. Instead, current approaches assume some form of local homogeneity for the material distribution, and/or uniformity in the stress state. In this manuscript, we describe a method that does not make these assumptions, but leads to a challenging inverse problem in OCE. Previously, we have demonstrated such a method in two dimensions, using the assumption of plain strain [12]. However, the direct extension of such a method to three dimensions is not tractable, due to the large number of optimization parameters. We demonstrate here, for the first time, its iterative solution in three dimensions. This solution is made computationally feasible by using adjoint-based algorithms [13], developing adaptive meshing techniques that focus computational efforts in regions with most significant spatial variations, and domain decomposition techniques to analyze different sub-regions in parallel and independently of each other. We note that Khalil et al. [14,15] proposed a method within which local homogeneity is not assumed; instead the problem is simplified by using prior information about the spatial distribution of mechanical properties. For most problems, especially in in vivo and in situ settings, this prior information is not available.
In this paper, we demonstrate our method in the context of compression OCE [16][17][18], in which the whole field of view being imaged undergoes quasi-static compression, and phasesensitive optical coherence tomography (OCT) is used to determine the resulting displacement field along the optical axis. We insert a transparent, compliant layer with known shear modulus between the sample and the compressor [18], and use this knowledge, along with the measured axial displacement in the layer and the sample, to recover the quantitative shear modulus map. The combined experimental and computational method to recover shear modulus distribution is applied to tissue-mimicking phantoms, and to ex vivo and in vivo human tissue samples. The tissue-mimicking phantom study highlights the ability of this method to recover accurate, absolute quantitative estimates and spatial variations of shear modulus. The ex vivo human tissue study, performed on a breast specimen containing invasive ductal carcinoma, demonstrates its capability to visualize cancerous tissue and its robustness in heterogeneous materials. Finally, the in vivo study, which is performed on human skin, points to the applicability of this method in resolving small, localized features. Overall, we demonstrate an important advance towards the accurate determination of the elasticity of soft tissue with micrometer-scale resolution over millimeter-scale regions.

Methods
In this section, we first describe the overall approach to generating quantitative volumetric images of mechanical properties and, thereafter, we describe each component of this approach in detail. In our approach, the whole sample volume to be imaged and the transparent, compliant layer are subjected to microscale, quasi-static compression [18]. Figure 1(a) shows the imaging setup sion embedde are pre-compr ator provides OCT is used t mm lateral fie (relative to th sensitivity (se are inputs to throughout th sion phantom the inverse pr shear modulu mesh element Fig. 1(c), are top panel of F method, whic [18]. and Sec accuracy affo the thresholde terparts for th as a reference  performed using a fiber-based Fourier-domain OCT system operating at a center wavelength of 835 nm and a 3-dB bandwidth of 50 nm. The axial and lateral resolution were measured to be 7.8 μm (in air) and 11 μm, respectively. For the skin data, OCE was performed using a fiber-based Fourier-domain OCT system operating at a center wavelength of 1300 nm and a 3-dB bandwidth of 170 nm. The axial and lateral resolution were measured to be 5.5 μm (in air) and 20 μm, respectively. The displacement sensitivity for both systems was measured to be <2.2 nm.
Prior to OCE acquisition, a compliant layer was placed on top of the sample. The imaging probe, consisting of a glass window attached to a piezoelectric actuator (PZT), was used to preload the layer and the sample, such that complete contact was achieved between the window, layer and sample surface within the OCT field of view. Further, the preload was used to ensure that the majority of the sample was under compressive strain. This preload was in the range 5-10% in the phantoms and 10-20% in the tissues.
After the initial preload, the PZT was used to provide sequential step-wise loading and unloading on the scale of a few micrometers. OCE acquisition was synchronized with the loading, and cross-sectional data (B-scans) were captured in each location in the field of view, sequentially, in the loaded and unloaded positions, assembling volumetric data. For the breast tumor, four acquisitions in multiple lateral locations were stitched together to obtain a larger field of view.

Displacement estimation
We denote an OCT A-scan as ex ( p ) ( ( ) ) ) ( I z a z z ιφ = , where ( ) I z is the complex-valued intensity associated with the pixel at depth z. Displacement data is obtained by acquiring two Ascans, one denoted 0 , where λ is the mean wavelength of the OCT system and n is the refractive index of the sample (assumed to be 1.4).

Algebraic method
The algebraic method for solving the inverse elasticity problem in compression OCE operates under the assumption of uniform and uniaxial compression. It was proposed and described in detail by Kennedy et al. [18], and we describe it briefly here. Under the assumption of uniaxial compression, shear modulus is given as: , where σ is the axial stress, and ε is the axial strain. ε is evaluated from the gradient of displacement using a weighted least mean square estimator [19]; σ is estimated from the strain of the compliant layer, and its wellcharacterized shear modulus at the particular preload. The stress is assumed to propagate uniformly in depth. The shear modulus becomes: is the tangent shear modulus and ( , ) c x y ε is the strain of the compliant layer. Here, (x, y, z) denote the lateral and axial coordinates.

Inverse problem solution
We partition the domain of interest, Ω, into the volume occupied by the compliant layer, c Ω , and the volume occupied by the sample, s Ω , and let ( , , ) u x y z denote the axial component of measured displacement. We model the compliant layer and the sample as incompressible elastic materials, with the stress tensor given by 2 pI where p is the pressure, I is the identity tensor, μ is the shear modulus, and ε is the strain tensor. Within the compliant layer, we measure the average axial strain through the thickness, ( , ) c x y ε , that is applied to ensure proper contact. This, in conjunction with the knowledge of the elastic behavior of the calibration material [18], allows us to calculate the tangent shear modulus, ( , ) c x y μ , within the compliant layer, at the particular preload.
The inverse problem we wish to solve is: given u and c μ , determine the spatial distribution of the shear modulus, ( , , ) x y z μ , in the tissue. This problem is supplemented by the constraint that the equations of equilibrium must be satisfied, that is 0 ∇⋅ = σ . We solve this problem as a minimization problem, where we seek to find the distribution μ that minimizes the functional [20]: In the equation above, the first term represents the data mismatch and the second term represents regularization. In the data mismatch term, u  is the axial component of the predicted displacement field, which is constrained to satisfy the equations of equilibrium for a linear, incompressible elastic solid, and T is a spatially varying weighting function that accounts for spatial dependence of noise in the measured displacement. The regularization term, which is of the total variation (TV) type, addresses the ill-posed nature of the inverse problem and penalizes large variations in the recovered shear modulus. It does so without regard to the steepness of the slope of these variations and is, therefore, ideally suited for recovering the shear modulus when working with sharp features; like the presence of a stiff tumor within soft glandular tissue. In this term, α is the regularization parameter, and c is another parameter that ensures that this term has continuous derivatives (with respect to μ) when 0 μ ∇ = .
We note that the shear modulus, μ, influences the data mismatch term indirectly via the predicted displacement, u  , since the spatial distribution of μ determines u  through the solution of the equations of equilibrium. Any change in μ leads to a change in u  , which changes the data mismatch term. The goal is to find the distribution of μ that minimizes the sum of the mismatch and the regularization terms. This is accomplished by representing both μ and u  using the standard finite element basis functions, and then using a gradient-based quasi-Newton algorithm to solve the minimization problem Eq. (1). The gradient vector, whose components represent the change in π due to change in the value of μ at a finite element vertex, is determined by solving the forward problem (equations of equilibrium) and an adjoint problem. At every iteration of the quasi-Newton algorithm, given the current estimate for the shear modulus distribution, the following tasks are performed: (1) A forward elasticity problem with the boundary conditions derived from the measured displacement is solved to generate the predicted displacement field; (2) An adjoint elasticity problem, which is driven by the difference between the measured and predicted displacement fields, is solved to determine the adjoint displacement field; and (3) Using the forward and adjoint displacement fields, the gradient vector is determined and supplied to the quasi-Newton method, which then returns with an improved estimate for the shear modulus field. These steps are repeated until convergence, which is obtained when the gradient vector is smaller than a specified tolerance, or when the change of the displacement mismatch over the last five iterations is smaller than a specified tolerance. We note that while this method is relatively new to OCE, it has been used quite extensively in ultrasound based elastography [13,20,27].
The algorithm described above is efficient; however, its cost can become prohibitive when solving problems in 3D on a fine grid. For example, if we were to apply this algorithm to the problem of a tissue-mimicking phantom with a single inclusion (as in Fig. 1), we would solve for 31 million optimization variables, which is very challenging. To address this challenge, we have incorporated two novel approaches. First, when solving the inverse problem, we make use of adaptive spatial resolution. In particular, we first solve the inverse problem on a coarse grid. Thereafter, we define a mesh size field that is inversely proportional to the absolute value of the gradient of the recovered shear modulus. This is used to generate a new finite element mesh that has finer resolution wherever the shear modulus varies sharply. We repeat this step several times and arrive at a final mesh whose resolution is optimal for the problem at hand. We terminate the process of mesh modification when we arrive at a mesh where the group of smallest elements has a mesh size that is comparable to the pixel size for the displacement field. The final mesh is fine where the gradients in μ are large, and coarse in other places. This procedure is shown in Fig. 1(c), where four successive steps of mesh refinement and the corresponding shear modulus distributions are shown for the tissue phantom considered in Section 3.1. As we move from one step to another, the mesh clusters near the interface between inclusion and the background where the finer spatial resolution is most needed. This process is also illustrated in Visualization 1, where the evolution of the shear modulus distribution and the underlying mesh is captured for the breast tissue example described in Section 3.2.
The approach described above works well for problems with a large field of view and relatively few significant features. In cases where there are many features within the field of view, this approach will lead to a fine resolution everywhere and, therefore, a very large computational problem. In order to tackle these types of problems, we combine this approach with a domain decomposition method. In this approach, the overall problem is simultaneously solved on N distinct meshes. In each mesh, an adaptive strategy is applied to a smaller subdomain of the total volume, and a fixed, coarse mesh is used in the rest of the domain. The fact that the adaptive strategy is applied to a single subdomain in each problem allows us to focus on the features contained within that subdomain, and the fact that the region outside this subdomain is represented through a coarse discretization, allows us to account for the effect of the remainder of the region without expending large computational resources. Another advantage of this method is that it is "trivially parallelizable"; that is, each of the N problems described above can be solved independently of each other on separate processors or computers, with no communication. In this manuscript, this domain decomposition approach was applied to the in-vivo imaging of the skin presented in Section 3.3 and is described in Fig. 5 for N = 3 subdomains.
The algorithms described above were implemented within our in-house finite element code written in Fortran 90. All problems were solved on multi-processor machine with eight 2.3 GHz CPU cores and a total memory of 16 GB. For a typical 3D problem with around 50,000 grid points and 400 iterations of the inverse algorithm, our code takes around 20 hours of wall clock time with the majority of this time spent solving the linear systems at each iteration.

Construction and characterization of phantom
Tissue-mimicking phantoms and the compliant layers were fabricated from two-compound room-temperature vulcanizing silicone. Two compounds were used: Elastosil P7676 and Elastosil RT601 (Wacker Chemie AG, Munich, Germany). The optical properties were controlled through the addition of titanium dioxide particles of known size distribution. The refractive index of the silicones is 1.4 and 2.5 for the titanium dioxide [21]. The mechanical properties were controlled by changing the mixing ratio of the two components of the silicone (Parts A and B), and further varied by diluting the silicone parts with polydimethylsiloxane (PDMS) oil (AK50, Wacker Chemie AG, Munich, Germany). PDMS oil was also used to lubricate the samples to minimize friction in the experiments. Each composition is presented as Compound A:B:Oil (e.g., P7676 2:1:0.3). The phantom mixture was mixed and sonicated to ensure uniform distribution of constituents and scatterers, and was oven-cured at 80°C.
The compliant layers were fabricated from P7676 1:1:0 with no added scattering. The layers were fabricated to a 500-μm thickness. The bulk of the four-inclusion phantom, described in Section 3.1, was fabricated from P7676 1:1:0, and the inclusions were fabricated from (in order of descending modulus): RT601 5:1:0, RT601 10:1:10, RT601 10:1:30 and P7676 2:1:0.3. 1 g/L of titanium dioxide was added to the bulk and 2.5 g/L was added to the inclusions to provide optical contrast. Each inclusion material was fabricated independently and then manually cut into an approximately rectangular shape, and incorporated into the bulk during its curing process. Each inclusion was designed to be 6 × , 3 × , 1.5 × and 0.5 × the bulk stiffness, with actual values shown in Fig. 2(a). Mechanical properties of each material were independently characterized using uniaxial compression testing (Instron, Norwood, MA, USA).

Preparation of breast tumor sample
Informed consent was obtained from patients and the study approved by the Human Research Ethics Committee of Royal Perth Hospital, Perth, Western Australia. The breast tumor sample was freshly excised and dissected for scanning, with approximate dimensions (x × y × z) of 2 × 2 × 1 cm 3 [18]. The sample was kept hydrated in saline until imaging, which occurred within 3 hours of excision. After imaging, samples were fixed in 10% neutral buffered formalin, embedded in paraffin, sectioned and stained with haematoxylin and eosin (H&E) following the standard histopathology protocols used at Royal Perth Hospital. Interpretation of histology was performed by an experienced pathologist.

Preparation of skin sample
Skin imaging was approved by The University of Western Australia Human Research Ethics Committee. Imaging was performed on a 26-year-old Caucasian male. Prior to imaging, the fingertip was cleaned and lubricated with PDMS oil. The compliant layer was also lubricated and placed on top of the fingertip. The finger was placed on a cushion to minimize discomfort and movement. The imaging probe was then brought into contact with the fingertip, preloading it until there were no observable movement artifacts and bulk motion in the displacement images. It was also preloaded to ensure that there was complete contact, and that all layers were undergoing strain due to compression. Figure 2 shows images of a phantom containing four inclusions spanning a range of shear moduli, including one that is less than that of the bulk matrix material. The raw OCT and shear modulus reconstructions are presented in both en face and B-scan images. The dashed line in the en face images indicates the position of the corresponding B-scan image. The primary purpose of this figure is to assess the performance of the reconstruction method on a sample whose mechanical properties have been measured through a standard compression test performed on a custom-built rig. We also compare the performance of the algebraic (Fig.  2(d)) and the iterative (Fig. 2(e)) methods. OCT images of the phantom (Fig. 2(a)) are presented in correspondence with the shear modulus images. Plots of the reconstructed shear modulus using the algebraic and iterative methods, along the white-dashed line on the en-face plane shown in Fig. 2(a), are shown in Fig. 2(b). Also, the percentage error in the average shear modulus for each inclusion and the surrounding material as a function of the modulus contrast (with respect to the background) is shown in Fig. 2(c) for the two methods. The mean and the standard deviation of the modulus within the inclusions and the background are reported in Table 1. These statistics were calculated by averaging the shear modulus over all mesh nodes within each region. The boundaries between the inclusion and the background were determined from the OCT volumes. The compression test was performed on separate samples that were created from the same mixture as the inclusions. Figure 2(c) demonstrates the accuracy of the iterative method in recovering the average shear modulus for all inclusions. The error is less than 10% for all inclusions, and the accura-cy is retained method grow tion of homog tion in accura as a fundamen method.

Ex vivo malignant breast tumor
Several studies have demonstrated that OCE can visualize tumor in excised human breast tissue [18,[23][24][25][26]. Eventually, OCE may provide intraoperative feedback to surgeons to assist in the detection of tumor, which may reduce the number of re-excision surgeries required to remove residual tumor missed during the initial surgery. In Fig. 3, we demonstrate our technique's ability to visualize tumor elasticity in human breast tissue acquired from a specimen containing invasive ductal carcinoma excised during a mastectomy procedure performed on a 60-year-old patient. Some of this data has previously been published based on the algebraic method [18]. The evolution of the shear modulus and the underlying mesh at various stages of the iterative inversion is shown in Visualization 1. The haematoxylin and eosin (H&E) histology image in Fig. 3 shows a large region of invasive tumor cells (T, stained purple) with local regions of benign ducts (D) in desmoplastic stroma. A region of stromal tissue (S), comprising a mixture of collagen, normal ductal structures, and microcystic structures, extends into the adipose tissue (A). The en face views in Figs. 3(c) and 3(d) show shear modulus images at a depth of 150 μm beyond the tissue-layer interface. The most striking difference between the images is the number of pixels with negative shear modulus value (marked by yellow) in the algebraic case, due largely to the problem of positive/tensile strain reported previously for OCE imaging of breast tissue [23]. In mechanically heterogeneous tissue, the algebraic method can result in strain in the direction opposite to the applied load, particularly at boundaries between tissue features. As the corresponding boundary stress is not considered by the algebraic method, this can lead to negative shear modulus values that are not physically meaningful. Furthermore, it is likely that these regions of non-physical shear modulus imply a reduced accuracy in the regions where the extracted modulus is physically plausible. In this case, depicted in Fig. 3(c), the pixels corresponding to regions of negative modulus are masked in yellow in the elastogram [18]. This problem is not present in the iterative method, where the shear modulus parameter is constrained to remain positive throughout the domain leading to a more accurate representation of the elasticity (Figs. 3(d)).
Although the elasticity is more faithfully estimated in adipose tissue using the iterative approach, the strongly fluctuating and low OCT signal typical of adipose tissue inevitably leads to an inaccurate result, as displacement accuracy and, therefore, elasticity accuracy are directly related to the OCT intensity. This could be the reason why adipose tissue and stroma show similar modulus values in the lower half of Fig. 3(d). In addition, in the region of tumor (T), a more uniform shear modulus is recovered using the iterative method, with the shear modulus in the range 60-65 kPa. By comparison, the corresponding region recovered using the algebraic method shows shear modulus values in the range 20-100 kPa. As the histology image demonstrates, there is relative uniformity in the structure of the tumor region, which is replicated in the shear modulus provided by the iterative method. It is worth noting that the uniform composition of the tumor here is distinct from many tumors previously scanned using compression OCE, where more commonly the intermingling of tumor with desmoplastic stroma creates heterogeneous strain as a hallmark of disease [26,27]. Also, in the histology image ( Fig. 3(b)), local regions of benign ducts in desmoplastic stroma are visible within the larger region of tumor. These regions are more clearly visible in shear modulus distributions obtained using the iterative solution by their lower modulus. The central duct descends past the tumor mass to the stromal region, whilst the others are circumscribed by the tumor. This is reflected in the iterative solution, suggesting that this approach can more faithfully reproduce features within mechanically heterogeneous tissue. images show t n ( Fig. 4(a)), fro ated keratinocy D), typically co n as spiraling s ermis; however dermis due to th show the shear layer (L)). It i more clearly, an Particularly, in so the modulus responding to F d e -E f = the ridges om top to ytes (coromprising structures r, there is he higher r modulus is evident nd further Fig. 4(i), s distribu- Fig. 4(c)) are readily id features are n The mean are: 34 kPa (S These average using the OC mean within e ods is challen strain, the eff such as inden gregate of a fe  -4(f)). The results from each refined partition in each inversion were fused to obtain a full, high-resolution reconstruction of shear modulus ( Fig. 4(g)).
Accurate mapping of mechanical properties in vivo in skin may find clinical application in dermatology in characterizing conditions such as cancer [30], pathological scarring [31] and scleroderma [32]. Further, it may find application in the cosmetics industry, where the mechanical properties of skin are important in aging, treatment and regeneration [33]. The first steps toward this with OCE have been recently demonstrated, showing the capacity to distinguish burn scars, moles and image other locations of the skin, however, all based on the algebraic method [34]. The iterative reconstruction method presented here will likely further improve the clinical diagnostic capacity of OCE.

Discussion
Optical coherence elastography has demonstrated great promise in mapping the mechanical properties of tissues on the hundreds of micrometer to millimeter scale, showing clinical potential, for example, in ophthalmology and oncology [8]. It has also started to emerge in cell mechanics, where the impetus is on understanding how cells interact with their mechanical microenvironment through physical forces [11]. The iterative inversion method presented in this paper will benefit all these applications of OCE leading to high-resolution and higheraccuracy shear modulus maps.
The iterative method to modulus reconstruction presented here makes no assumptions regarding homogeneity of stress or elastic properties. As a result, it produces modulus distributions that are more accurate when compared with the algebraic method, especially for samples that are inherently heterogeneous. As shown in the Results section, in the case of a phantom with multiple inclusions, it leads to more accurate average modulus values, especially when the contrast between the inclusion and the background is large. Also, consistent with how the sample was manufactured, it leads to a more uniform distribution within each inclusion. In the case of the ex vivo imaging of a malignant breast cancer specimen, the iterative method avoids the appearance of regions of negative modulus, and produces modulus images that are qualitatively more consistent with histology. Finally, in the case of in vivo imaging of skin, the iterative method better resolves the three layers of the skin (stratum corneum, epidermis and dermis), and the spiral structure of individual sweat ducts.
The improved accuracy of the iterative method comes at the cost of solving a complex inverse problem which is computationally intensive. These computational costs are mitigated by the use of an adjoint solver, an adaptive grid resolution algorithm, and a domain decomposition algorithm. The computational gains from this relatively simple approach are substantial. For example, in the multiple inclusion problem, using this approach, the number of unknown parameters in the inverse problem reduces from ~28.9 million to ~42 thousand, a reduction by a factor of 680 that makes the solution of this problem feasible in about 6 hours on an off-the-shelf desktop computer. Similar savings are observed for other problems reported in this paper.
It is interesting to consider the spatial resolution of the modulus images and the least count error for the modulus estimate produced by the iterative method. However, these quantities are not easily estimated because the modulus is determined from the displacement field by solving a complex nonlinear inverse problem on an adaptive unstructured grid, and the displacement estimate is itself obtained from phase differences between collocated pixels in successive OCT B-scans. One approach to estimating the spatial resolution would be to determine the point-spread function (PSF) for each of these steps and then convolve them in order to determine a cumulative PSF for the entire problem. However, for the complex, nonlinear inverse problem we are solving, it is not possible to determine the PSF. Thus, in order to gain some understanding of the spatial resolution and the least count error of the method, we must rely on some other approach. In particular, we make use of the validation results obtained for the multiple inclusion problem described in Fig. 1.
We assume that the interface between the inclusion and the background is infinitesimally thin, and that this discontinuity is smeared out in the modulus image. By measuring the risedistance, which is defined as the distance taken for the modulus value to rise from 10% to 90% of its value, we can obtain an estimate of the spatial resolution of the method. This distance was calculated for the inclusion in Fig. 1, and is equal to 62 μm. In order to estimate the least count error, we assume that the modulus distribution within the inclusion is homogeneous. This is a reasonable assumption since the modulus of the inclusions was controlled by the mixing ratio of the PDMS. Further, the addition of scatterers did not show a change in overall mechanical properties when assessed by uniaxial compression tests. Thus, any variation reported in the modulus value can be ascribed to an error in the reconstruction. The root mean square (RMS) of this error is, therefore, an estimate of the least count error. The RMS error value for the inclusion is ~4 kPa.
It is interesting to consider the effect of overall noise, as well as it spatial variation, on the results of the iterative method, especially since we expect the noise in the displacement measurement to vary spatially in inverse proportion to the OCT signal intensity. Any increase in the overall level of noise leads to an increase in the regularization parameter through the Lcurve procedure, which in turn leads to smoother reconstructions with lower contrast due to the effect of the TV regularization.
The spatial variation in noise is accommodated through the spatially varying weighting function T in Eq. (1). This function determines the relative importance of data mismatch and regularization terms at any point in the domain. In regions with large noise this function is set to a small value thereby effectively increasing the importance of regularization. When viewed from a Bayesian inference perspective, with a Gaussian model for noise, this function is the inverse of the covariance of noise in the measured data. To understand what effect this might have on a reconstruction, let's consider two regions: one where the noise estimate is small (T is large) and another where it is large (T is small). In the region with small noise (large T) the regularization term will have a lesser effect and the modulus reconstruction will have significant spatial variations. In the region where noise is large (small T), the regularization term will be more important and the solution will be smoother. In the context of the examples considered in this manuscript, where we have varied T in inverse proportion to OCT intensity, this implies that regions with lower OCT intensity will have relatively smoother reconstruction when compared with regions with higher OCT intensity.
While the method described in this paper does not require the assumption of homogeneity that is invoked in the algebraic method, it retains some other assumptions, namely linear elasticity, and isotropy. For the phantom and breast samples considered in this paper, isotropy is likely a reasonable assumption. However, it is also likely that given the amount of precompression, these samples are not in the linear-elastic regime. Thus, the modulus we recover can be interpreted as a tangent modulus about the pre-compression strain, which can be characterized by the absolute deformation state of the compliant layer. A natural extension of the iterative method discussed here would be to include non-linear elastic behavior, anisotropy, and perhaps even poroelastic behavior in the inverse problem. We note that the iterative method provides a framework to include all of these effects [20,[35][36][37]. However, to successfully do so will require collecting additional independent displacement data that may be obtained by altering the spatial and/or temporal characteristics of the excitation, and enriching the forward solver within the iterative method with the appropriate mechanical models.

Conclusion
We have demonstrated, for the first time, a tractable method for the measurement of shear modulus in three dimensions on length scales highly relevant to biological and disease processes on the scale between cells and whole tissues. This approach is free from simplifying assumptions that have been invoked in other methods and are routinely violated in heterogeneous samples. Undoubtedly, there is still much to learn in the application of these methods to a wide variety of biological constructs and tissue types. However, the results demonstrated here suggest that our approach will find many applications in the future, in research and with prospects for translation to applications from tissue regeneration to clinical medical imaging.

Funding
Australian Research Council; the Cancer Council, Western Australia; and the Department of Health, Western Australia. PM is supported by a Royal Society University Research Fellowship.