Analysis of spatial resolution in phase-sensitive compression optical coherence elastography

Optical coherence elastography (OCE) is emerging as a method to image the mechanical properties of tissue on the microscale. However, the spatial resolution, a main advantage of OCE, has not been investigated and is not trivial to evaluate. To address this, we present a framework to analyze resolution in phase-sensitive compression OCE that incorporates the three main determinants of resolution: mechanical deformation of the sample, detection of this deformation using optical coherence tomography (OCT), and signal processing to estimate local axial strain. We demonstrate for the first time, through close correspondence between experiment and simulation of structured phantoms, that resolution in compression OCE is both spatially varying and sample dependent, which we link to the discrepancies between the model of elasticity and the mechanical deformation of the sample. We demonstrate that resolution is dependent on factors such as feature size and mechanical contrast. We believe that the analysis of image formation provided by our framework can expedite the development of compression OCE. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Optical coherence elastography (OCE) is emerging as a higher resolution alternative to both ultrasound elastography and magnetic resonance elastography in a range of applications including in oncology [1][2][3], cardiology [4,5], ophthalmology [6,7], and dermatology [8,9]. Although the higher spatial resolution of OCE has been purported to be a main advantage over the more established elastography techniques, it has yet to be clearly defined or measured. As a result, it is challenging to evaluate and compare the performance of OCE systems, and to identify the most suitable applications or, indeed, which variant to use in a given application.
In medical imaging more generally, spatial resolution is often defined as the smallest distance that two objects can be brought together and still be distinguishable as separate objects, and is an important metric to assess image quality. In optical coherence tomography (OCT) and other optical microscopy techniques, resolution is treated as a system parameter, i.e., independent of the sample, and can be estimated from the point-spread function (PSF) of the imaging system [10,11]. Following the precedent set by OCT, until now, OCE resolution has been loosely defined as a combination of system parameters, namely, the resolution of the underlying OCT system, and the signal processing used to map displacement or velocity to a mechanical property, such as elasticity [12][13][14][15]. However, in OCE, the resolution of mechanical properties is intrinsically linked to the mechanical deformation of the sample. Whilst this has been suggested in the literature, it has yet to be investigated or demonstrated.

Background
OCE elastogram generation can be described in three stages; the deformation of the sample under mechanical loading, the measurement of this deformation using OCT, and the signal processing used to relate measured deformation to elasticity. In this section, we describe each of these stages in the context of phase-sensitive compression OCE.

The contribution of mechanical deformation to resolution
Typically, the goal of elastography is to recover the mechanical properties of the sample by applying mechanical loading and observing the resulting deformation. In compression OCE, the sample is placed under a static preload (referred to as the unloaded state) and a microscale compressive load is then applied (referred to as the loaded state). The change in position of each point in the sample is described by the displacement field, u. The sample in both the unloaded and loaded states is taken to be in a state of static equilibrium and, additionally, the deformation between the two states is assumed to be much smaller than the size of the sample [30]. As a result, the equations of equilibrium reduce to: where σ is the infinitesimal Cauchy stress tensor field, which describes the internal stress field present in the sample [31]. The resulting deformation is described by the corresponding infinitesimal strain tensor, ε, whose components, ε ij , are derived from the partial derivatives of u [32]: Commonly, the sample is assumed to be a linear elastic solid [33], whose mechanical properties are represented by a fourth-order elasticity tensor, C, that relates the stress and strain tensors:
where μ and λ are the Lamé parameters and δ ij is the Kronecker delta, which is equivalent to 1 if i = j and 0 otherwise. In phase-sensitive compression OCE, only the axial component of the displacement field, u z , is measured [34]. Given an observed u z , to solve for the mechanical properties at each point in the sample, one would need to solve for σ. However, this requires challenging numerical methods, thus, simpler models of σ are typically employed, i.e., that σ is uniaxial and uniform throughout the sample [12][13][14][15]. Under these assumptions, elasticity is specified by Young's modulus (E), as a combination of the Lamé parameters, or the ratio of uniaxial stress, σ zz , over uniaxial strain, ε zz , (3 2 ) .
zz zz E σ μ λ μ λ μ ε That is, the elastogram of uniaxial strain is assumed to be inversely proportional to elasticity [12][13][14][15]. These assumptions have enabled elastography to be relatively easily and robustly performed with straightforward mechanical loading. However, sample characteristics such as

The contribution of OCT to resolution
In our implementation of compression OCE, axial displacement is measured using phasesensitive OCT. OCT images of the sample are modeled as a convolution of the sample scatterers with the complex OCT PSF [10]. The OCT resolution, and hence the resolution of the displacement measurement, are defined by the width of the PSF envelope. For this study, we consider 2D elastograms presented in the x-z plane where the envelope of the OCT PSF is defined as [11]: where w z and w x are the 1/e 2 radii of the OCT intensity in z and x respectively; and I 0 is a scaling factor. By convention, OCT resolution is typically given in terms of the full-width-athalf-maximum (FWHM) of the intensity PSF where the axial and lateral resolution is 2ln 2 times w z and w x , respectively.

The contribution of least-squares regression to resolution
For the elastograms presented in this study, least-squares linear regression (LSR) is used to estimate the axial strain as the gradient of axial displacement with depth, over a sliding window of length Δz, typically 15 μm -100 μm [12][13][14] to alleviate the impact of system noise [34]. As all measurements of axial displacement are acquired with equal spacing, ordinary least-squares linear regression is equivalent to convolution with a Savitzky-Golay (SG) kernel [35]. A 2n + 1 point, first order, SG differentiation kernel is given by [35]: where β = 3/n(n + 1)(2n + 1) is the smoothing coefficient corresponding to fitting a first degree polynomial (i.e., a line) [36] and h is the spacing of the voxels in the axial dimension. Convolution of SG(z) with a unit step displacement, a(z), gives the PSF corresponding to LSR: For a fit length given by Δz = h(2n + 1), the FWHM of the PSF for LSR is: In this study, the system resolution in compression OCE is modeled by the FWHM of the result of convolving PSF OCT (x,z), Eq. (6), with PSF LSR (z), Eq. (8). Feature resolution is modeled as the convolution of the system resolution with the effective local PSF from mechanical deformation determined from FEA, as detailed in Section 3.4.

Simulation
To simulate OCE, we combine a finite-element model of mechanical deformation using the Abaqus simulation software package (Dassault Systèmes, Providence, USA, v6.14), with a model of the strain, given phantom, sho linear-elastic were fixed on coefficient of lower limit of and experime with 50 μm q the feature sid simulations, t square inclusi changing the h As the ph sample defor displacement case is simula solving for th variation in m feature, the m phantom, the the choice of measure the s simulation is impact of def in Fig. 2

Tissue-s
To validate s silicone elast were 3 mm and AK50 Sil using titanium OCT system g in Eq. (8). T own in Fig. 2 [40]. The mechanical properties were controlled by selecting different elastomers and varying the ratio of the catalyst, curing agent, and silicone oil. The stress-strain relationship of each silicone was characterized using a custom-built uniaxial compression testing apparatus. The 2D plane-strain model used in the simulations effectively assumes that the sample is infinitely long in the out-of-plane (y) dimension [39] and to best match this assumption, the inclusions were made longer (10 mm long rectangular prisms) in the y-dimension, than in the x-and z-dimensions. The inclusions were cut using a custom made silicone slicing tool, to improve the accuracy of the cross-section geometry. Each phantom was made in four steps to ensure the inclusion was placed 600 μm from the top surface. Firstly, a base layer of silicone was poured into a glass petri dish, using real-time OCT in the x-z plane to aid in achieving the correct height, which, for each inclusion size, was determined as the total phantom height (3 mm) minus the sum of 600 μm and the inclusion height. Secondly, once the base had cured, the inclusion was placed on top. Thirdly, silicone was poured over the inclusion and base, guided by real-time OCT, to a height of 600 μm above the inclusion. Finally, once cured, the phantom was removed from the petri dish using a 15 mm diameter biopsy punch. Five different phantoms were made with varying inclusion size and mechanical contrast, as shown in Table 1.

Phase-sensitive compression OCE system
OCE measurements were performed using a fiber-based spectral-domain OCT system (Telesto 320, Thorlabs Inc., USA). The light source is a superluminescent diode with a mean wavelength of 1300 nm and a spectral bandwidth of 170 nm. The measured OCT axial resolution in air is 4.8 μm (FWHM of irradiance). The objective lens (OL in Fig. 3(a)) (LSM03, Thorlabs) has a measured lateral resolution in air of 7.2 µm (FWHM of irradiance). AK50 silicone oil was applied to lubricate the phantom-imaging window interface. A preload of 15% bulk strain was applied to each phantom using a translation stage to ensure uniform contact between the rigid plate (RP in Fig. 3(a)), the phantom (P in Fig. 3(a)), and the imaging window (IW in Fig. 3(a)). The ring actuator (RA in Fig. 3(a)) (Piezomechanik GmbH, Germany) has an aperture of 65 mm and a maximum stroke of 10 µm. A 75 mm diameter imaging window (Edmund Optics Inc., USA), fixed to the ring actuator, transfers the compressive load from the actuator to the phantom. The system was operated in commonpath where the imaging window itself, a partial reflector, acted as the OCT reference reflection. Scans were acquired by taking 2,000 A-scans per B-scan, and 2,000 B-scans per C-scan over a 5 mm × 5 mm region with a lateral sampling density of 2.5 μm per voxel. The ring actuator was driven in a quasi-static regime by a 12.5 Hz square wave, collinearly with the imaging beam and synchronized with the acquisition of OCT B-scans, such that alternate B-scans are acquired at different compression levels. Local axial displacement, u z , in Fig. 3(b), is calculated from the change in phase, Δφ(z), between B-scans acquired at the same lateral location [41]: The WLS weight to me ratio (SNR) [3 least-squares operation, in in Section 2 experiment ar dimension.
In experim contact betwe applied aroun this microscal

Quantify
We measure f feature bound of which is de illustrated in elastograms o the feature bo 4(b), respectiv an error funct which is fit to scaling factor he mean wave roximately 1.4 cement using 34]. 3

Effect of
In Fig. 6, we axial and late phantoms of 250 μm Figs. 6(a)-6(c stiffness (at 1 and lateral fe denoted by th In simula to 1000 μm and lateral fe respectively. These results d ut instead varie ne, like many s e entire sample mation. The lat usion and fricti se boundaries, t in axial strai ion boundaries of the bulk is l xial inclusion b nts decrease to e resolution in t 5 As the fea feature, restri at the feature with the scale strain to be lin of the feature above the fea increased, the will become i region, limitin be seen by the the lateral edg phantom, res in Fig. 6(h). m six indepen nd the axial Figs. 6(g) and and 6(h) show to feature size, nd 100 μm to 4 6 Fig. 7(g) iffness increas extent, increa f the inclusion increases, the radient in axia he latter effe r this feature g ng mechanical f system resol illustrate the r on. We analyz resolution and nd lateral resolu resolution was Phantoms 4 and to 250 µm. It i ession will sm the approach for Phantom 1 50 μm to 20 50 μm to 250 μ were generated plane (y) dimen axis of the mean ± stand ed from simula om FEA alone limit to feature 1, 4 and 5, res 8 Fig. 8, as the axial system resolution is improved, the feature resolution approaches the limit imposed by mechanical deformation. For Phantom 1, the feature resolution was reached to within 10% of the deformation limit with a system resolution of 36 μm (Δz = 50 μm). For Phantoms 4 and 5, the feature resolution was reached to within 10% of the deformation limit with a system resolution of 107 μm (Δz = 150 μm). This result demonstrates that prior knowledge, or reasonable assumptions about the sample geometry, can enable the optimal system resolution to be designed for a given application.

Discussion
We have presented the first study of resolution in OCE, by combining a finite-element model of mechanical deformation, with a model of the OCT system and a model of the signal processing, both based on linear systems theory, validated by close correspondence with experiment. This study illustrates the limitations imposed on resolution by sample mechanics and the elastogram formation process. Previously, Kennedy et al. analyzed the impact of mechanical deformation on contrast in OCE, and showed that sample-specific factors such as geometry and mechanical contrast limit elastogram accuracy and sensitivity [45]. Their study employed FEA to model deformation of elastic features in the absence of system noise to investigate the theoretical limits of elastogram contrast. Separately, Chin et al. combined FEA with a linear systems model of OCT to demonstrate that, in addition to the commonly considered effects of optical noise, the coupling between mechanical deformation and its optical detection impacts the precision of phase-sensitive compression OCE [46]. A promising avenue for future work is to integrate these previous studies into a similar framework presented in this study to develop a more complete model of image formation in OCE.
The results presented in Section 4 demonstrate that due to the limit imposed by mechanical deformation, feature resolution in compression OCE varies along feature boundaries and is dependent on feature size and mechanical contrast. Improving the OCT and signal processing resolution beyond this deformation limited state will not improve the feature resolution, but will instead unnecessarily trade off other image quality metrics such as contrast, sensitivity, imaging depth and field-of-view. The proportionality between feature resolution and size is predominately due to the scale invariant nature of deformation [44] and the incompressibility of the sample, leading to non-uniform stress that is not accounted for in the mechanical model. Similarly, the dependence of feature resolution on mechanical contrast is also attributed to incompressibility. Due to the assumption of uniform stress and the incompressibility of many soft tissues [43], similar results could be expected in tissue.
A square inclusion geometry was used in this study to readily measure a step response in both axial and lateral directions. From the theory presented in Section 2 and the results presented in Section 4, a different feature shape, e.g., a cylindrical or spherical inclusion, will likely exhibit a different intra-sample resolution distribution. As shown in Section 4.1, the incompressibility of the sample and adhesion of the soft bulk to the stiff inclusion leads to gradients in axial strain at the inclusion boundaries, degrading feature resolution. At the corners of the rectangular inclusions, the deformation of the bulk is less restricted, effectively improving axial and lateral resolution in these locations. In the case of cylindrical or spherical inclusions, which effectively have a circular cross-sectional geometry, the axial deformation of the bulk at the inclusion boundaries will be restricted in a different way, likely resulting in a different distribution of feature resolution. The focus of this study was to demonstrate that the sample has a pronounced effect on the resolution of OCE, and to establish a framework that incorporates the deformation of the sample with that of the OCT and signal processing to analyze resolution in compression OCE; however, an analysis of the effect of feature shape on resolution represents an avenue for future research.
The assumptions in the mechanical model used to relate deformation to elasticity in OCE can directly affect feature resolution, and importantly, suggests that a more complex mechanical model that makes fewer assumptions about sample deformation can improve feature resolution. This can be achieved by using numerical solutions to the inverse elasticity problem as demonstrated by Dong et al. [29]. Their work employs an iterative approach to solving Eqs. (1) and (4) for the mechanical properties of the sample and, importantly, removes the assumption of uniform and uniaxial stress. Preliminary results from elastograms presented in their study demonstrate a homogenization of resolution across the field-of-view that is achieved by reducing the link between feature properties and resolution. However, the improvement to feature resolution through this approach is yet to be quantified. Further removing assumptions of sample homogeneity, linearity and isotropy can likely lead to improved feature resolution in OCE, where our framework could be used to quantify these improvements.
A potential application of OCE is in tissue engineering and cell mechanics [47]. Contemporary developments in this area seek to mimic the mechanical complexity of the extracellular environment, moving towards complex, multi-parameter, 3D biomaterials [48]. OCE is well-poised to provide a non-invasive technique that can characterize biomaterials in 3D. Previous studies have used compression OCE systems with an OCT axial and lateral resolution of below 2 μm to study cellular scale structures in animal models and excised tissue [25][26][27]. For example, features with a lateral length of 20 μm in excised mouse aorta were observed with a measured equivalent Gaussian FWHM lateral feature resolution of approximately 10 μm to 15 μm [27], which is consistent with our framework. However, the diameter of many animal cells lies below 10 μm [49], and to study cells with a feature resolution below this limit will likely require solving more complex mechanical models, or indeed employing techniques that localize deformation [19][20][21][22]. The limit on achievable resolution may result in compression OCE being most useful in studying the mechanical properties over the entirety of the cell and the interactions between the cell and its environment over relatively large fields-of-view. Our framework can help to address this question.
In the simulations presented in this study, we consider the OCT and least-squares regression resolution independently, however, in experiment, the resolution of the OCT system determines speckle size [50]. Due to the impact of speckle on the accuracy of displacement measurements [30] and therefore strain accuracy [34], the fit length of axial displacement used in least-squares regression is effectively dependent on the OCT resolution. This dependence between OCT and least-squares regression resolution was not considered in our study, but could be incorporated into our framework by combining previous work by Chin et al. [46]. From the results in this study, the degradation to feature resolution from least-squares regression was typically much greater than OCT. Our framework suggests that improving displacement accuracy in phase-sensitive OCE can lead to significant improvements in feature resolution by removing the need for larger fit lengths in least-squares regression.
This study has focused on phase-sensitive compression OCE, however, our framework could be extended to other elastography techniques. For example, quantitative microelastography (QME), an extension of compression OCE, uses a compliant silicone layer placed between the sample and the imaging window to map the 2D uniaxial stress applied at the sample surface [51]. The stress is assumed to be uniform with depth, which, along with the volume of local axial strain, provides an estimate of Young's modulus throughout the sample. However, as QME uses the distribution of local axial strain to determine the distribution of Young's modulus, the spatial resolution of QME is tethered to the resolution of compression OCE. Consequently, feature resolution in QME will also likely be spatially varying and sample dependent where our framework could be applied to assess feature resolution. Further, alternative OCE loading methods, such as shear wave and surface acoustic wave OCE [6,17,52], typically rely on the same assumption of homogeneity to access sample mechanical properties. It is likely that mechanical heterogeneity will lead to similar discrepancies between feature and system resolution. A study comparing all such methods with our framework would be informative and the feature resolution could help determine the most efficient approach for various applications. Similarly, our framework can be extended to quantify improvements in feature resolution brought about by high resolution OCT systems [25][26][27], improved signal processing [28,29] and extended to shear wave OCE and surface acoustic wave OCE or indeed to ultrasound elastography [53] and magnetic resonance elastography [54].

Conclusion
We have presented the first framework to study resolution in phase-sensitive compression OCE that incorporates the mechanical deformation of the sample in response to a load, the measurement of the resulting displacement in the sample using OCT, and the signal processing used to estimate local strain. This framework enables us to analyze the impact of these factors on both system and feature resolution. We have demonstrated, through simulation and experiment, that mechanical deformation is typically the dominant determinant of resolution in compression OCE, and imposes a limit to resolution. Due to this limit imposed by mechanical deformation, the ability to resolve features in an underlying distribution of mechanical properties varies along feature boundaries and is dependent on feature size, mechanical contrast and location.

Funding
Australian Research Council; the Cancer Council Western Australia; OncoRes Medical; William and Marlene Schrader Trust of the University of Western Australia scholarship.