Analysis of image formation in optical coherence elastography using a multiphysics approach

Image formation in optical coherence elastography (OCE) results from a combination of two processes: the mechanical deformation imparted to the sample and the detection of the resulting displacement using optical coherence tomography (OCT). We present a multiphysics model of these processes, validated by simulating strain elastograms acquired using phasesensitive compression OCE, and demonstrating close correspondence with experimental results. Using the model, we present evidence that the approximation commonly used to infer sample displacement in phase-sensitive OCE is invalidated for smaller deformations than has been previously considered, significantly affecting the measurement precision, as quantified by the displacement sensitivity and the elastogram signal-to-noise ratio. We show how the precision of OCE is affected not only by OCT shot-noise, as is usually considered, but additionally by phase decorrelation due to the sample deformation. This multiphysics model provides a general framework that could be used to compare and contrast different OCE techniques. ©2014 Optical Society of America OCIS codes: (000.3860) Mathematical methods in physics; (000.4430) Numerical approximation and analysis; (030.6140) Speckle; (110.2990) Image formation theory; (110.4500) Optical coherence tomography. References and links 1. B. F. Kennedy, K. M. Kennedy, and D. D. Sampson, “A review of optical coherence elastography: fundamentals, techniques and prospects,” IEEE J. Sel. Top. Quantum Electron. 20(2), 1–17 (2014). 2. J. M. Schmitt, “OCT elastography: imaging microscopic deformation and strain of tissue,” Opt. Express 3(6), 199–211 (1998). 3. B. F. Kennedy, S. H. Koh, R. A. McLaughlin, K. M. Kennedy, P. R. T. Munro, and D. D. Sampson, “Strain estimation in phase-sensitive optical coherence elastography,” Biomed. Opt. Express 3(8), 1865–1879 (2012). 4. A. Nahas, M. Bauer, S. Roux, and A. C. Boccara, “3D static elastography at the micrometer scale using full field OCT,” Biomed. Opt. Express 4(10), 2138–2149 (2013). 5. M. Razani, A. Mariampillai, C. Sun, T. W. H. Luk, V. X. D. Yang, and M. C. Kolios, “Feasibility of optical coherence elastography measurements of shear wave propagation in homogeneous tissue equivalent phantoms,” Biomed. Opt. Express 3(5), 972–980 (2012). 6. S. Song, Z. Huang, T.-M. Nguyen, E. Y. Wong, B. Arnal, M. O’Donnell, and R. K. Wang, “Shear modulus imaging by direct visualization of propagating shear waves with phase-sensitive optical coherence tomography,” J. Biomed. Opt. 18(12), 121509 (2013). 7. C. Li, G. Guan, X. Cheng, Z. Huang, and R. K. Wang, “Quantitative elastography provided by surface acoustic waves measured by phase-sensitive optical coherence tomography,” Opt. Lett. 37(4), 722–724 (2012). 8. J. Li, S. Wang, R. K. Manapuram, M. Singh, F. M. Menodiado, S. Aglyamov, S. Emelianov, M. D. Twa, and K. V. Larin, “Dynamic optical coherence tomography measurements of elastic wave propagation in tissuemimicking phantoms and mouse cornea in vivo,” J. Biomed. Opt. 18(12), 121503 (2013). 9. S. G. Adie, X. Liang, B. F. Kennedy, R. John, D. D. Sampson, and S. A. Boppart, “Spectroscopic optical coherence elastography,” Opt. Express 18(25), 25519–25534 (2010). #212910 $15.00 USD Received 27 May 2014; revised 17 Jul 2014; accepted 23 Jul 2014; published 1 Aug 2014 (C) 2014 OSA 1 September 2014 | Vol. 5, No. 9 | DOI:10.1364/BOE.5.002913 | BIOMEDICAL OPTICS EXPRESS 2913 10. W. Qi, R. Li, T. Ma, J. Li, K. Kirk Shung, Q. Zhou, and Z. Chen, “Resonant acoustic radiation force optical coherence elastography,” Appl. Phys. Lett. 103(10), 103704 (2013). 11. V. Crecea, A. Ahmad, and S. A. Boppart, “Magnetomotive optical coherence elastography for microrheology of biological tissues,” J. Biomed. Opt. 18(12), 121504 (2013). 12. R. K. Wang, Z. Ma, and S. J. Kirkpatrick, “Tissue Doppler optical coherence elastography for real time strain rate and strain mapping of soft tissue,” Appl. Phys. Lett. 89(14), 144103 (2006). 13. R. K. Wang, S. Kirkpatrick, and M. Hinds, “Phase-sensitive optical coherence elastography for mapping tissue microstrains in real time,” Appl. Phys. Lett. 90(16), 164105 (2007). 14. C. Sun, B. Standish, B. Vuong, X.-Y. Wen, and V. Yang, “Digital image correlation-based optical coherence elastography,” J. Biomed. Opt. 18(12), 121515 (2013). 15. J. Fu, F. Pierron, and P. D. Ruiz, “Elastic stiffness characterization using three-dimensional full-field deformation obtained with optical coherence tomography and digital volume correlation,” J. Biomed. Opt. 18(12), 121512 (2013). 16. S. G. Adie, B. F. Kennedy, J. J. Armstrong, S. A. Alexandrov, and D. D. Sampson, “Audio frequency in vivo optical coherence elastography,” Phys. Med. Biol. 54(10), 3129–3139 (2009). 17. B. F. Kennedy, T. R. Hillman, R. A. McLaughlin, B. C. Quirk, and D. D. Sampson, “In vivo dynamic optical coherence elastography using a ring actuator,” Opt. Express 17(24), 21762–21772 (2009). 18. B. F. Kennedy, M. Wojtkowski, M. Szkulmowski, K. M. Kennedy, K. Karnowski, and D. D. Sampson, “Improved measurement of vibration amplitude in dynamic optical coherence elastography,” Biomed. Opt. Express 3(12), 3138–3152 (2012). 19. P. D. Ruiz, J. M. Huntley, and J. M. Coupland, “Depth-resolved imaging and displacement measurement techniques viewed as linear filtering operations,” Exp. Mech. 51(4), 453–465 (2011). 20. S. Song, Z. Huang, and R. K. Wang, “Tracking mechanical wave propagation within tissue using phase-sensitive optical coherence tomography: motion artifact and its compensation,” J. Biomed. Opt. 18(12), 121505 (2013). 21. A. S. Khalil, R. C. Chan, A. H. Chau, B. E. Bouma, and M. R. Mofrad, “Tissue elasticity estimation with optical coherence elastography: toward mechanical characterization of in vivo soft tissue,” Ann. Biomed. Eng. 33(11), 1631–1639 (2005). 22. A. S. Khalil, B. E. Bouma, and M. R. Mofrad, “A combined FEM/genetic algorithm for vascular soft tissue elasticity estimation,” Cardiovasc. Eng. 6(3), 93–102 (2006). 23. K. M. Kennedy, C. Ford, B. F. Kennedy, M. B. Bush, and D. D. Sampson, “Analysis of mechanical contrast in optical coherence elastography,” J. Biomed. Opt. 18(12), 121508 (2013). 24. G. Lamouche, B. F. Kennedy, K. M. Kennedy, C.-E. Bisaillon, A. Curatolo, G. Campbell, V. Pazos, and D. D. Sampson, “Review of tissue simulating phantoms with controllable optical, mechanical and structural properties for use in optical coherence tomography,” Biomed. Opt. Express 3(6), 1381–1398 (2012). 25. M. M. Doyley, P. M. Meaney, and J. C. Bamber, “Evaluation of an iterative reconstruction method for quantitative elastography,” Phys. Med. Biol. 45(6), 1521–1540 (2000). 26. R. C. Chan, A. H. Chau, W. C. Karl, S. Nadkarni, A. S. Khalil, N. Iftimia, M. Shishkov, G. J. Tearney, M. R. Kaazempur-Mofrad, and B. E. Bouma, “OCT-based arterial elastography: robust estimation exploiting tissue biomechanics,” Opt. Express 12(19), 4558–4572 (2004). 27. V. Crecea, A. L. Oldenburg, X. Liang, T. S. Ralston, and S. A. Boppart, “Magnetomotive nanoparticle transducers for optical rheology of viscoelastic materials,” Opt. Express 17(25), 23114–23122 (2009). 28. V. Y. Zaitsev, L. A. Matveev, A. L. Matveyev, G. V. Gelikonov, and V. M. Gelikonov, “Elastographic mapping in optical coherence tomography using an unconventional approach based on correlation stability,” J. Biomed. Opt. 19(2), 021107 (2014). 29. V. Y. Zaitsev, L. A. Matveev, G. V. Gelikonov, A. L. Matveyev, and V. M. Gelikonov, “A correlation-stability approach to elasticity mapping in optical coherence tomography,” Laser Phys. Lett. 10(6), 065601 (2013). 30. S. Makita, F. Jaillon, I. Jahan, and Y. Yasuno, “Noise statistics of phase-resolved optical coherence tomography imaging: single-and dual-beam-scan Doppler optical coherence tomography,” Opt. Express 22(4), 4830–4848 (2014). 31. T. A. Krouskop, T. M. Wheeler, F. Kallel, B. S. Garra, and T. Hall, “Elastic moduli of breast and prostate tissues under compression,” Ultrason. Imaging 20(4), 260–274 (1998). 32. B. F. Kennedy, R. A. McLaughlin, K. M. Kennedy, L. Chin, A. Curatolo, A. Tien, B. Latham, C. M. Saunders, and D. D. Sampson, “Optical coherence micro-elastography: mechanical-contrast imaging of tissue microstructure,” Biomed. Opt. Express 5(7), 2113–2124 (2014). 33. K. M. Kennedy, S. Es’haghian, L. Chin, R. A. McLaughlin, D. D. Sampson, and B. F. Kennedy, “Optical palpation: optical coherence tomography-based tactile imaging using a compliant sensor,” Opt. Lett. 39(10), 3014–3017 (2014). 34. J. M. Schmitt and A. Knüttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A 14(6), 1231–1242 (1997). 35. A. Curatolo, B. F. Kennedy, D. D. Sampson, and T. R. Hillman, “Speckle in Optical Coherence Tomography,” in Advanced Biophotonics: Tissue Optical Sectioning, V. V. Tuchin, and R. K. Wang, eds. (Taylor & Francis, 2013), pp. 211–277. 36. C.-E. Bisaillon, G. Lamouche, R. Maciejko, M. Dufour, and J.-P. Monchalin, “Deformable and durable phantoms with controlled density of scatterers,” Phys. Med. Biol. 53(13), N237–N247 (2008). #212910 $15.00 USD Received 27 May 2014; revised 17 Jul 2014; accepted 23 Jul 2014; published 1 Aug 2014 (C) 2014 OSA 1 September 2014 | Vol. 5, No. 9 | DOI:10.1364/BOE.5.002913 | BIOMEDICAL OPTICS EXPRESS 2914 37. J. M. Schmitt, A. Knüttel, and R. F. Bonner, “Measurement of optical properties of biological tissues by lowcoherence reflectometry,” Appl. Opt. 32(30), 6032–6042 (1993). 38. F. J. van der Meer, D. J. Faber, D. M. Baraznji Sassoon, M. C. Aalders, G. Pasterkamp, and T. G. van Leeuwen, “Localized measurement of optical attenuation coefficients of atherosclerotic plaque constituents by quantitative optical coherence tomography,” IEEE Trans. Med. Imaging 24(10), 1369–1376 (2005). 39. C. Xu, J. M. Schmitt, S. G. Carlier, and R. Virmani, “Characterization of atherosclerosis plaques by measuring both backscattering and attenuation coefficients in optical coherence tomography,” J. Biomed. Op


Introduction
Optical coherence elastography (OCE) is a technique in which an image (elastogram) is formed of a mechanical property of a sample. In OCE, a mechanical load is applied to a sample, and the resulting deformation is detected using optical coherence tomography (OCT) [1]. There are many forms of OCE, differentiated by the method of mechanical loading, such as compressive [2][3][4], shear wave [5,6], surface wave [7,8], frequency-swept [9,10] and localized loading using magnetic nanoparticles [11]. Similarly, there are several methods for measuring the resulting sample deformation, including phase-sensitive detection [12,13], speckle tracking [2,4,14,15], and Doppler spectrum analysis [16][17][18]. Regardless of the implementation, inherent to OCE is the interaction between two processes: the mechanical deformation of the sample and its detection using OCT. A theoretical framework for describing both processes is vital in understanding elastogram formation, and in establishing the performance of OCE on a variety of samples and for different system parameters. Several groups have previously analyzed aspects of the elastogram formation process, including recent studies focused on the detection of sample displacement [3,19,20] and earlier studies examining the mechanical deformation of samples [21][22][23]. These studies have largely considered the processes of deformation and detection as independent.
In this paper, we present evidence that the coupling between mechanical deformation and its detection by OCT significantly affects the measurement precision of OCE (i.e., the sensitivity to which the sample deformation is detected), to an extent not previously considered. We do this through a multiphysics simulation of OCE, which combines a finite element model of mechanical deformation, capable of simulating complex geometries, and a linear systems model of displacement detection by OCT, incorporating attenuation and shot-noise-limited optical detection. Multiphysics simulation has the advantage that it can model a wider range of sample properties, system parameters and loading conditions than is generally represented with OCE phantom studies [24]. Previous studies have proposed multiphysics models for simulating speckle-tracking based techniques in both ultrasound elastography [25] and OCE [26]. In this study, we demonstrate a multiphysics simulation approach to phase-sensitive OCE, which is the most prevalent signal processing method currently used [1]. We validate this model by simulating elastograms acquired using a specific OCE technique, namely phase-sensitive compression OCE [3,13], and demonstrate close correspondence with experimental results. We then use this model to analyze the measurement precision of phasesensitive OCE. We show that the displacement sensitivity, and hence the precision, is affected not only by optical noise, as is generally considered [3,10,12,13,23,27], but additionally by "phase decorrelation noise", which results from a violation of the assumption, when measuring displacement from OCT scans, that sample deformation preserves phase correlation. Phase decorrelation noise varies throughout an OCE elastogram as a function of the local deformation of the sample. Previous studies on speckle-tracking methods in OCE have shown that there is a threshold level of decorrelation beyond which sample displacement cannot be reliably estimated [15,26,28,29]. Similarly, a recent study on Doppler flow imaging analyzed the adverse effect of decorrelation on Doppler phase sensitivity [30]. However, to the best of our knowledge, the effect of decorrelation has not previously been considered in phase-sensitive OCE. We demonstrate that existing techniques of measuring displacement sensitivity in phase-sensitive OCE [3,10,12,13,23, 27] overestimate system performance as the amount of sample deformation is increased. Additionally, we demonstrate how this variation in measurement precision affects elastogram quality by using our model to analyze how strain signal-to-noise ratio (SNR) varies with sample deformation. Whilst we focus our analysis on phase-sensitive compression OCE, similar relations between displacement sensitivity and decorrelation may hold true in other forms of OCE. In principle, this framework is extendable to other forms of mechanical loading and OCT-based detection methods, providing a tool for comparing and contrasting the variants of OCE.

Phase-sensitive optical coherence elastography
In phase-sensitive OCE, the axial displacement within a sample, d(x,y,z) in response to a load, is calculated from the change in the OCT phase, Δφ(x,y,z), between scans of the loaded and unloaded sample [3,12,13], i.e., where λ 0 is the mean free-space wavelength of the OCT system, and n is the sample refractive index at location (x,y,z). In compression OCE, the axial displacement can be used to calculate the local axial strain (i.e., the strain defined over a finite range) [3,12,13] and in shear wave and surface acoustic wave OCE, it can be used to calculate the phase velocity of the propagating wave [5][6][7][8].

Phase-sensitive compression OCE
This study focuses on phase-sensitive compression OCE. This technique requires only two OCT scans: one of the sample unloaded, or under a static preload, and another of the sample under an additional compressive load described by an applied stress tensor, σ [1]. The resulting sample deformation can be quantified by the strain tensor, ε [1]. For strains < 0.1, it has been reported that tissues such as breast and prostate, and elastomers such as silicone, can be described by a linear elastic model [24,31,32]. Additionally, considering the sample to be isotropic, if the compressive load is uniform and uniaxial at the sample surface, then the uniaxial stress, σ z , and strain, ε z , of the sample are related through a single constant, Young's modulus, E = σ z /ε z , which is a common measure of material stiffness. In general, the Young's modulus will vary throughout the sample, and a map of Young's modulus can be calculated by measuring the local strains and stresses within the sample. The local strain can be determined from the spatial derivative of the measured displacement at each point within the sample. It is not straightforward, however, to measure the local stress distributed throughout the material, although our group has recently proposed a technique to measure local stress at the tissue surface [33]. Compression elastograms are typically, thus, a map of the strain distribution throughout the sample in response to the load. Phase-sensitive detection generally measures only the axial displacement, so the local, axial strain, ε z (x, y, z), in phase-sensitive compression OCE is calculated as where Δd(x,y,z) is the change in axial displacement over the depth range Δz. In this study, the phase difference, Δφ(x,y,z), is unwrapped prior to calculating the axial displacement to increase the measurable dynamic range [32]. The axial strain is then calculated using weighted-leastsquares (WLS) linear regression over a sliding window on the axial displacement, with the weights being the OCT signal-to-noise ratio of irradiance (OCT SNR) associated with each displacement measurement [3].

Metrics of elastogram quality and precision
We consider the precision of OCE as quantified by three metrics: displacement sensitivity, strain sensitivity, and strain SNR. In phase-sensitive OCE, the minimum measurable displacement is determined by the phase difference sensitivity of the OCT system, s Δφ . Assuming shot-noiselimited detection, and OCT SNR >> 1, s Δφ ≈(SNR OCT ) −1/2 [3]. The displacement sensitivity due to optical noise, s dO , is then The typical method of calculating the displacement sensitivity of an OCE system is to first measure the phase difference sensitivity (the standard deviation of the phase difference) from OCT scans of a stationary, unloaded reflector, then apply Eq. (3) [3, 10, 12, 13, 23, 27]. However, because of phase decorrelation noise, this does not quantify the true displacement sensitivity in OCE. Phase decorrelation noise is caused by both strain-induced and translationinduced decorrelation. Strain-induced decorrelation, φ dε , arises from mechanical deformation, or strain, of the sample, which causes the sub-resolution scatterers in a particular region to move closer or further apart, which in turn causes decorrelation of the OCT speckle pattern, and its phase, introducing additional errors into the displacement measurement. Translation-induced decorrelation, φ dt , arises from shifts in the mean location of scatterers in a particular region due to the loading. The true displacement sensitivity, s d , is then a function, f d (), of three terms, In this study, we measure the displacement sensitivity of a sample under varying loads. The method of measuring displacement sensitivity described above works well for point measurement techniques, such as Doppler OCT; however, in OCE, a collection of displacement measurements acquired from different spatial locations is required to estimate the mechanical parameter of interest [1]. For example, compression OCE requires displacement measurements from a number of axial points to calculate the local strain [1,3]. To take this into account, and to more accurately define displacement sensitivity in the context of compression OCE, in this study we use a measure for displacement sensitivity by considering multiple points within an OCT scan which are all undergoing the same displacement but have a range of OCT SNR values. In a homogeneous sample subject to uniform, uniaxial compression, all lateral points at a given depth undergo the same axial displacement. Thus, we calculate displacement sensitivity as the standard deviation of the measured displacement over a lateral extent at a given depth in the loaded sample. As the relationship between displacement sensitivity and OCT SNR is nonlinear, the displacement sensitivity of a collection of points with varying OCT SNR, with a mean of SNR µ , is less than the displacement sensitivity of a single point measurement with OCT SNR equal to SNR µ . The displacement varies with depth, so taking measurements at different depths allows us to analyze the effects of translation-induced decorrelation on the displacement sensitivity. Similarly, all points in a homogeneous sample will experience the same axial strain, so taking measurements under different loads allows us to analyze the effects of strain-induced decorrelation on the displacement sensitivity. In addition, the strain sensitivity is measured as the standard deviation of the strain over a mechanically homogeneous region. The strain SNR is then given by the ratio of the mean strain, µ ε , over the strain sensitivity, s ε , expressed using a log scale as [3,23] ( )

Multiphysics model of optical coherence elastography
In this study, the precision of phase-sensitive OCE is analyzed using a multiphysics simulation model. The model we present comprises four components: simulation of the optical image formation using a linear systems model of OCT, simulation of the sample mechanical deformation using the finite element method (FEM), combining the mechanical deformation with the detection by OCT, and generation of the OCE elastogram.

OCT image formation
The model of OCT image formation presented here is based on previous linear systems theory formulations [34,35]. The sample is represented as a collection of discrete scattering potentials; this has been shown to be a good approximation to many tissue-mimicking phantoms for OCT and OCE, which are typically fabricated by embedding discrete scattering particles in a medium [24,[34][35][36]. The OCT image is approximated by convolving a complex point-spread-function (PSF) with the map of scattering potentials. Finally, noise is added under the assumption of shotnoise-limited detection. The scattering potentials (represented by dimensionless impulses) which comprise the sample are distributed randomly over a discrete grid with sub-wavelength spacing in the axial and lateral dimensions given by δz i and δx i , respectively. The relative density of scattering potentials within different regions of the model was chosen to match the relative concentration of scatterers within corresponding regions in the sample, e.g., the physical phantom. The magnitudes of the scattering potentials are scaled to represent the backscattering coefficient of each scattering particle, which is constant in this study, as we assumed a constant particle size and type, and low numerical aperture (NA). The magnitudes are then further scaled according to the Beer-Lambert law of attenuation, which is commonly used in OCT to model the attenuation of samples [37-41]. The backscattered OCT irradiance, I(z), as a function of depth, z, is proportional to the reflectance, R(z) = ρexp(−2µ t z), where ρ is the backscattering coefficient, and µ t is the attenuation. The amplitude of the backscattered OCT electric field per unit area, A(z), is then proportional to the square root of the reflectance, proportional change in reflected field amplitude, Using Eq. (6) and µ t (x,y,z), a map of the attenuation throughout the sample, we can build an attenuation-scaled map of the reflected field amplitude, RM(x,y,z), given by Multiplying the map of scattering potentials by RM(x,y,z) then gives an attenuation-scaled map of scattering potentials, H(x,y,z).
In this optical simulation, the backscattering and attenuation coefficients of the scattering potentials can be varied arbitrarily. Physically, the two will be related based on the type, size, and concentration of scatterers in the sample. If the type and size of scattering particles is the same throughout the sample, then both the OCT backscattered irradiance and the attenuation coefficient are proportional to the particle concentration [36].
Details of the OCT complex PSF are given elsewhere [34, 35], but are summarized here for completeness. Supposing a Gaussian source spectrum, Gaussian beam illumination and collection, a low numerical aperture (NA < 0.1), negligible sample-induced aberration and negligible contribution from multiple scattering, a 3-D PSF, Ψ(r,z), can be defined based on radial coordinate r (where r 2 = x 2 + y 2 ) and axial coordinate z, with amplitude A(r,z) and phase Φ(r,z) given, respectively, by Here, is the coherence length of the source with free-space mean wavelength λ 0 , 3-dB bandwidth Δλ and mean sample refractive index n, assuming no dispersion over the source bandwidth; , with a s and R being the 1/e radii of the irradiance of the Gaussian beam at the focus, and before entering the objective lens, respectively; z 0 = kR 2 is the effective Rayleigh range of the unfocused beam in the medium; is the mean wavenumber in the medium; f = f 0 n is the focal length in the sample of the objective lens with free-space focal length f 0 ; L = f − z is the physical distance from the lens; and . The simulated complex OCT scan is evaluated by convolving Ψ(r,z) with H(x,y,z), and sampling the result at the axial and lateral intervals specified by the OCT system, δz j and δx j , respectively. This is implemented by directly evaluating the sum, S(x j ,y j ,z j ), at every voxel (x j ,y j ,z j ) in the output, i.e., x y = + , and (x i ,y i ,z i ) is the location of the i-th scattering potential. This implementation avoids having to specify a new discrete sampling grid for the scattering potentials when they are displaced under FEM-computed deformation (Section 3.3), thus, allowing resolution of the scatterer locations to arbitrary precision, whilst restraining the computational memory requirements.
Finally, the OCT system is taken to be operating in the shot-noise limit, with shot-noise simulated by adding isotropic, independent, complex Gaussian white noise (with noise power, σ 2 ) [43, 44]. Figure 1 shows an example of the simulated OCT image formation process demonstrated using a structured, tissue-mimicking phantom [45], with comparisons to experimental results. Figures 1(i) and 1(k) demonstrate the close correspondence between simulated and measured OCT images.

Mechanical deformation
Mechanical deformation of the sample under an applied load is modeled using the finite element method, a numerical method for computing approximate solutions to boundary value problems by subdividing the problem into a finite number of discrete, homogeneous elements [46]. The coupled equilibrium equations are then solved at each element, given a set of model parameters [46]. A FEM model is constructed by: 1. Defining the geometry of the sample; 2. Assigning material properties to each of the deformable regions of the sample; 3. Subdividing the model into discrete elements -this is referred to as "meshing" the model, and the resulting elements as the "FEM mesh"; 4. Applying boundary conditions, such as surface friction; and 5. Defining a known load or displacement introduced to the sample, corresponding, in this case, to the load applied to the sample during OCE imaging.
The solution of the FEM model then provides the displacements, strains and stresses experienced by each of the elements in the mesh. Specifically, the displacements are evaluated at the vertices of the mesh elements, and the strains and stresses are evaluated at the integration points which lie within the mesh elements [46]. Figure 2 shows an example of a FEM simulation of compressive loading applied to a sample comprising a stiff rectangular inclusion embedded, at a slight angle to the surface, in a softer bulk surrounds. The size of the FEM mesh, and the size of the displacements, has been exaggerated for visibility.

Combining optical and mechanical models
Simulating OCE requires application of the mechanical deformations from the solution of the FEM model to the map of optical scattering potentials to generate an OCT simulation of the sample under load. We accomplish this using linear interpolation, with barycentric coordinates, of the locations of the scattering potentials. In particular, the displacements calculated from the FEM model, which are known only at the vertices of the FEM mesh, are interpolated throughout the computational space to all locations where scatterers reside.
The interpolation method uses, in general, a different mesh to that used by the FEM model, allowing the FEM model to be solved using an element shape which is optimal for the sample geometry. Once the FEM solution is obtained, Delaunay triangulation [47] is used to obtain a second mesh, upon which the interpolation is based, which uses the same vertices as the FEM mesh but simplex elements (triangles in 2-D, tetrahedra in 3-D). Delaunay triangulation is a common triangulation method with optimized implementations in many programming environments. Barycentric interpolation is then used to linearly interpolate the calculated displacements from the vertices of the FEM mesh to any required location.
For ease of notation, we consider the case of a 2-D field of scattering potentials, however, the same methodology applies in 3-D. The location of a scattering potential, r s = (x, z) T , can be expressed in terms of a linear combination of the three vertices (r 1 , r 2 , r 3 ) of the triangle which surrounds the scattering potential, r s = λ 1 r 1 + λ 2 r 2 + λ 3 r 3 . The real-valued coefficients λ 1 , λ 2 , λ 3 are called the barycentric coordinates of r s with respect to (r 1 , r 2 , r 3 ), and fulfill the constraint that Given a general function f() that transforms r 1 , r 2 , r 3 to f(r 1 ), f(r 2 ), f(r 3 ), the barycentric interpolation of f(r s ) is then given by f(r s ) = λ 1 f(r 1 ) + λ 2 f(r 2 ) + λ 3 f(r 3 ), where the barycentric coordinates λ 1 , λ 2 , λ 3 are the same as those computed from Eq. (11). Let f() be the function which takes r i , the location of the i-th mesh vertex in the unloaded sample, and returns f(r i ), the location of the same vertex in the loaded sample; f(r i ) = r i + u i , where u i is the FEM computed displacement of the i-th vertex from the unloaded to the loaded sample. Applying this transformation to the location, r sj , of the j-th scattering potential in the sample then gives the location, f(r sj ), of the scattering potential in the loaded sample. Figure 3 shows a schematic of this process. Applying Eq. (10) to the new locations of the scattering potentials, it is then possible to generate a simulated OCT scan of the sample under the applied load. A flowchart illustrating the simulation process is shown in Fig. 4. The inputs and simulation parameters are summarized in Table 1. The outputs are the simulated OCT scans of the sample before and after loading, and the simulated OCE elastogram.

Phantom fabrication and characterization
To validate the OCE simulation, a cylindrical tissue-mimicking phantom (diameter ∼2.5 cm, thickness ∼1 mm) with controlled optical and mechanical properties was fabricated using silicone elastomers (Wacker, Germany) [23, 24, 36]. The phantom comprised a stiff inclusion in softer bulk surrounds. The optical scattering properties were controlled using known concentrations of titanium dioxide (TiO 2 ) particles (mean diameter 1 µm) evenly mixed into the silicone (at concentrations of 2.5 mg/mL in the inclusion and 0.8 mg/mL in the bulk). The mechanical properties were controlled by selecting different elastomers and varying the ratio of the catalyst (A) to the curing agent (B). The inclusion was made using Elastosil RT601, in the ratio 5:1 (A:B), and the bulk using Elastosil P7676, at 2:1 (A:B). The mechanical properties of the silicones were characterized using a standard compression testing apparatus (Instron, Norwood, Massachusetts). In addition, to enable analysis of the effects of mechanical deformation on the precision of OCE measurements, a homogeneous phantom was simulated using properties similar to those of the soft bulk silicone in the inclusion phantom.

Compression OCE system and optical parameters
OCE measurements were performed using a fiber-based Fourier-domain OCT system [3,23,32]. The light source is a superluminescent diode with a mean wavelength of λ 0 = 835 nm, and a 3-dB spectral bandwidth of Δλ = 50 nm. The objective lens in the sample arm has a focal length of f 0 = 35.8 mm in air, and the 1/e radius of the intensity of the Gaussian beam before the lens was measured to be R = 793 μm. The measured free-space axial/lateral resolution (full width-at-half maximum irradiance) is 8.5 µm/11 µm. The group refractive index of the samples is approximately n = 1.42. The interference spectrum for each A-scan is captured over 1,300 pixels of a CMOS line-scan camera, resulting in a free-space axial imaging range of 3.13 mm, a freespace axial sampling density of 4.8 µm/pixel, and an effective imaging range and sampling density in the sample of 2.20 mm and δz j = 3.4 µm/pixel, respectively. The lateral sampling density was δx j = 4 µm/pixel. Scans were acquired with a lateral range of 4 mm, cropped to 2.5 mm.
The sample attenuation coefficient, μ t , was measured from experimental scans of the silicone phantoms, and σ 2 , the variance (power) of the optical shot-noise, was chosen such that the resulting average OCT SNR matched the experimental scans in selected corresponding areas. The relative concentrations of the scattering potentials in the model were chosen to match the ratios of TiO 2 particles used in the different parts of the phantoms, and the initial positions of the scattering potentials were discretized to a grid spacing of δz i = δx i = 0.5 μm. For the inclusion phantom, the attenuation was measured to be ∼3.69 mm −1 in the inclusion and ∼1.18 mm −1 in the bulk. The average OCT SNR was 18.6 dB in an area 325 µm × 50 µm (x × z) at the top of the phantom. The relative density of scattering potentials in the inclusion:bulk was 3.125:1. In the homogeneous phantom, the attenuation coefficient was 1.5 mm −1 . The noise level was chosen to produce an average OCT SNR of 15 dB at the top of the homogeneous phantom.
The sample arm contained an imaging window rigidly fixed to a piezo-electric ring actuator, enabling imaging and mechanical loading of the sample from the same side [17, 23]. A preload, generating 9-22% bulk strain, was applied to each sample using a rigid brass plate, of larger surface area than the sample, to ensure uniform contact between the brass plate, the sample, and the imaging window [23]. The system operated in a common-path configuration, and the imaging window itself, a partial reflector, acted as the OCT reference arm mirror.
The piezo-actuator was driven by a 5 Hz square wave, slow enough that the sample could be considered to be under quasi-static load. This was synchronized to the OCT B-scan acquisition rate of 10 Hz, ensuring consecutive B-scans were acquired in the loaded/unloaded state [3]. Axial displacements in the sample were calculated from Eq. (1) using the unwrapped phase difference between consecutive (loaded minus unloaded) OCT B-scans. A schematic of our phase-sensitive compression OCE setup, and example displacement and strain A-scans are shown in Fig. 5. The zero-phase reference coincides with the position of the imaging window (labeled "reference reflector" in Fig. 5(a)) in this common-path setup. The phase difference, Δφ, and hence the relative displacement between the reference reflector and the sample, d, are both zero at the imaging window, and maximal at the rigid plate and the measured displacement and strain values are negative under compressive loading. The local strain was calculated using weighted-leastsquares linear regression over a fitting range of 100 μm [3].

Mechanical parameters
The FEM models were constructed in the Abaqus simulation software package (Dassault Systèmes, Providence, USA, v6.12). The geometry of the inclusion phantom was determined from structural OCT scans of the phantom. The elasticity of the samples was modeled using a linear, isotropic, elastic model. Using the stress/strain curves obtained from the bulk compression testing of the constituent phantom materials, a value for Young's modulus was calculated from a tangent fitted to the curve about a quiescent point specified by the bulk strain applied during the preload stage of imaging [23]. The inclusion phantom was subject to a preload displacement of ∼270 µm from an initial thickness of ∼930 µm, corresponding to an initial pre-strain of 22%. This gives a Young's modulus of 830 kPa in the inclusion, and 33.3 kPa in the bulk, comparable to inclusion phantoms examined in previous studies [23]. Both the inclusion and the bulk are modeled as being nearly incompressible, with Poisson's ratio close to 0. 5 [23], and the displacement introduced to the surface of the phantom by the piezo-actuator is 0.93 µm. The homogeneous phantom was simulated with an initial pre-strain of 9% and a thickness under preload of ∼1060 µm, giving a Young's modulus of 19.2 kPa, comparable to homogeneous phantoms examined in previous studies [23]. The Poisson's ratio is also taken to be close to 0.5, and the displacement introduced to the surface of the phantom was simulated for 23 discrete steps between 0.21 µm and 15 µm.
Since the phantoms are relatively wide compared with the imaging field-of-view (FOV), we approximated the 3-D sample deformation by a 2-D plane strain model, which set the displacement and strain in the y direction to zero [46]. Similarly, lateral symmetry conditions were defined at the left and right edges of the FEM model, which set lateral (x) displacement at the edges of the OCT FOV to zero. The silicone phantoms were relatively sticky (likely due to a portion of the silicone remaining uncured), and not lubricated during OCT scanning, so friction at the top and bottom surfaces of the FEM simulation was taken to be infinite (no lateral motion at these surfaces) during loading. The average spacing between vertices of the FEM mesh was set to ∼5 µm. Initial tests showed that this gave comparable numbers, with much shorter processing time, than finer mesh spacings. Figure 6 shows the simulation of phase-sensitive compression OCE compared to experimental scans of the silicone inclusion phantom. Overall, the results show good agreement between the experiment and simulation. Variations in the axial strain can be seen surrounding the inclusion (Figs. 6(e) and 6(f)) due to non-uniform stresses caused by the inclusion geometry [23]. Figure 7 shows five regions used for numerical comparisons of the mean OCT SNR, mean displacement (µ d ), displacement sensitivity (s d ), mean strain (µ ε ), strain sensitivity (s ε ), and strain SNR (SNR ε ) between the experimental and simulated scans. The results of these comparisons are shown in Table 2. OCT and strain measurements were averaged over each entire region; displacement and displacement sensitivity were evaluated along a line at the bottom of each region. The numbers are comparable between the experiment and simulation, although the displacement and strain sensitivity are less in the experimental scans in all cases. These results are discussed in more detail in Section 6.

Displacement sensitivity
Displacement sensitivity was estimated using simulations of the homogeneous phantom detailed in Section 4. If calculated in the usual manner, from multiple measurements of a single point with high OCT SNR in a stationary, unloaded sample, the displacement sensitivity is 2.2 nm from a point with an OCT SNR of 27 dB. This is comparable to previous experimental measurements of 1.2 nm at an OCT SNR of 50 dB [23], which, although calculated at a higher OCT SNR, is less than would be expected from Eq. (3) due to phase noise caused by environmental conditions and galvanometer mirror noise. As detailed in Section 2.2, in this study, we calculate the displacement sensitivity as the standard deviation of the displacement over a 500 µm lateral extent, as this better reflects that OCE elastograms generally require multiple phase or displacement measurements to estimate strain [1]. Figure 8 shows plots of the displacement sensitivity, calculated as defined above, as a function of strain and depth in the sample. For this homogeneous sample, the axial strain is constant throughout the sample for a given load. The axial displacement is likewise constant at any given depth within the sample; increasing depth corresponds to an increase in the relative displacement between the loaded and unloaded sample with respect to the reference reflector (see Fig. 8(b)). The average OCT SNR at the top of the phantom was 15 dB (min ∼0 dB, max ∼27 dB). Even in the absence of noise (blue curves), the displacement sensitivity degrades, due to phase decorrelation, with both increasing strain ( Fig. 8(a)) and depth ( Fig. 8(b)) into the sample. In both plots, the displacement sensitivity decreases dramatically as the combined phase decorrelation (strain-induced and translation-induced) reaches a threshold.
The black lines show the displacement sensitivity calculated using the method described in Section 2.2, from the simulation of a stationary, unloaded sample; the sensitivity was calculated to be 14.8 nm at a mean OCT SNR of 15 dB. This value provides an upper bound on the plots of displacement sensitivity in the presence of both optical noise and phase decorrelation (red curves). There is a strain threshold below which the displacement sensitivity is limited by optical noise (solid red curve in Fig. 8(a)); however, decorrelation with depth always causes the displacement sensitivity to decrease, regardless of the initial displacement sensitivity. Figure 9 shows plots of the strain SNR calculated over 500 µm × 50 µm (x × z) regions in the simulated homogeneous phantom. With no optical noise (blue curves), strain SNR is constant with mean strain at a particular depth, until a threshold is reached at which the strain SNR decreases sharply ( Fig. 9(a)). When optical noise is included, increasing strain increases strain SNR until the same threshold is reached. Similar to the displacement sensitivity, even without including optical noise, the strain SNR decreases with depth into the sample, until a point at which it drops abruptly ( Fig. 9(b)). The higher the applied strain, the greater is the axial displacement at a particular depth, and the closer is the point of strain estimation failure to the sample surface.

Discussion
This study has demonstrated the first multiphysics simulation of phase-sensitive OCE. Related models have previously been proposed for simulating speckle-tracking based methods in OCE.   In this paper, we used our simulation to assess the effect of sample deformation on the precision of measured displacement and strain. A related study by our group examined the accuracy of OCE strain elastograms in representing the elasticity distribution of samples [23]. That study showed how uneven stress distributions surrounding a stiff inclusion can lead to "wings" of high strain around the inclusion, which matches what we observe in Figs. 6(e) and 6(f). However, the previous study did not consider the interaction between sample deformation and detected optical signal. In this paper, we have demonstrated that this interaction sets bounds on measurable displacement and strain in phase-sensitive OCE. A major implication, not considered in the previous study, is that the precision of OCE strain elastograms, and hence the precision to which we can infer the sample elasticity distribution, is heavily dependent on the value of the applied load.
It should be noted that the phase decorrelation noise described in this study is not noise in the conventional sense of being random with each measurement, but rather, like speckle, it is the deterministic result of a particular arrangement of sub-resolution scatterers. Measuring displacement with phase-sensitive OCE is only valid under the assumption that the speckle pattern in the scans used to determine the phase shift is fully correlated, analogous to the frozen speckle model described previously by Duncan and Kirkpatrick [49]. Phase decorrelation noise is then the degradation in precision that results when this assumption is no longer valid. Consider the case of a single scattering particle at the focus of the OCT beam, as it moves a small distance from z = 0 to z = δz. From Eq. (9), the phase shift resulting from this motion is ( ) ( ) In Section 5.1, we demonstrated the ability of our multiphysics simulation framework to model elastograms generated using phase-sensitive compression OCE. Scans of the inclusion phantom shown in Fig. 6 demonstrate good agreement between the experimental data and the simulation; however, the speckle patterns clearly differ between the two OCT B-scans in Figs. 6(a) and 6(b). This is due to the infeasibility of matching the exact locations of the sub-resolution scattering particles in the experiment and the simulation [35]. The current optical simulation also does not take into account effects such as the confocal function, or wave-front distortion in the turbid medium, which may be present in the experimental scans. Nevertheless, provided that the statistics of fully developed speckle patterns match, as is the case here, the differences in the specific speckle realizations do not affect the validity of the multiphysics simulation framework in modeling elastograms. In Table 2, the displacement sensitivity and strain SNR are slightly inferior in the experimental data, acquired from scans of a 3-D phantom, compared to the 2-D simulation. This is likely caused by sample motion out-of-the-plane (along the y-axis) in the experimental scan that could result in additional phase decorrelation. The largest discrepancies between the experiment and simulation are in the regions closest to the boundaries (top and bottom) of the phantom. This is likely due to the effects of friction [23], which is assumed to be infinite in the simulation, but in the experimental scan is likely to be finite; the exact value is, however, unknown.
Section 5.2 demonstrates the effects of phase decorrelation on the displacement sensitivity of phase-sensitive OCE. In the absence of optical noise, the displacement sensitivity degrades approximately linearly with the amount of strain in the sample ( Fig. 8(a)), although at lower strains the displacement sensitivity is limited by the optical noise rather than by strain-induced decorrelation. However, from Fig. 8(a), it is clear that this best-case sensitivity degrades rapidly when the sample is subject to loading, i.e., the displacement error is increasing with strain. Without phase unwrapping of the phase difference, Δφ, the maximum detectable strain for a 1 mm thick sample is << 1 mε, which, combined with lower OCT SNR in real samples, is likely the reason previous studies on OCE have not encountered this issue. The translation-induced decorrelation is likely to be a more significant factor in practice. As Fig. 8(b) demonstrates, the displacement sensitivity degrades approximately exponentially with depth into the sample, which, in our system, corresponds to increasing relative displacement between the reference reflector and the sample. Translation-induced decorrelation degrades the displacement sensitivity regardless of the OCT SNR; beyond a depth of ∼800 µm into the sample, the displacement sensitivity can have degraded by a factor in the range of 2 to 5. Section 5.3 shows the corresponding effects of phase decorrelation on the strain SNR (Fig. 9). Even without optical noise, strain-induced decorrelation limits the maximum strain SNR to approximately 40 dB for this system, which occurs close to the zero-phase reference at the top of the phantom. Without noise, the ideal scenario is reached with low levels of applied strain, as these achieve the best displacement and, hence, strain sensitivity, whilst retaining the maximum achievable strain SNR at the surface. With noise, there is a tradeoff between increasing strain SNR with higher applied load, and achieving better displacement and strain sensitivity with lower applied load. This suggests that for optimum strain SNR in phase-sensitive compression OCE, the applied load should be below the strain threshold, but otherwise as large as possible.
The results shown in Fig. 8 and Fig. 9 were acquired from simulations of homogeneous samples in which the strain resulting from uniform compressive loading is uniform across the sample. In more complex geometries, such as the inclusion geometry shown in Fig. 6, or in tissue, the strain can vary greatly throughout the sample. In particular, the strain in areas surrounding a stiff inclusion, such as in a phantom or around a region of fibrosis in tissue, can be much higher than the strain in the surrounding bulk. Nevertheless, the general conclusions to be drawn from the results shown in Fig. 8 and Fig. 9 still hold, as can be seen in Table 2. Regions 1 and 2 in Table 2 have comparable OCT SNR, but Region 1 has a higher mean strain, just past the threshold of ∼2 mε, and therefore a moderately lower strain sensitivity, but still higher strain SNR than Region 2, as would be expected; similarly for Regions 2 and 3.
A common method to improve OCT SNR is to average multiple scans. Improvement by averaging assumes that the acquired data consists of a noise component, which is isotropic and independent with each measurement, and a signal component, which is constant over all the averaged points. These conditions hold in the case of OCT shot-noise, but not for decorrelation noise, which, for the same loading conditions, does not change between OCT measurements. Thus, averaging of the OCT scans will improve the shot-noise contribution to displacement and strain sensitivity, bringing the red curves of Fig. 8 and Fig. 9 closer to the blue ones, but we expect it not to have any effect on the decorrelation noise contribution. A possible method to minimize decorrelation noise could be to sample the displacement at multiple points during loading, such that the displacement, and hence relative strain, between any two sample points is minimized. However, this would require many more OCT acquisitions and, thus, slow the OCE scan speed.
This study has focused on phase-sensitive compression OCE; however, the multiphysics simulation framework we have presented here is readily extendable to other forms of mechanical loading and optical detection. For example, needle OCE [50, 51] could be simulated by restricting the simulation output to a single A-scan. Dynamic loading methods such as shear wave [5,6], surface wave [7,8], or frequency-swept loading [9,10], could be modeled by running the simulation for each time step. Similarly, the simulation framework is extendable to modeling deformation in 3-D.

Conclusion
We have presented a multiphysics model for elastogram formation in OCE, which incorporates the mechanical deformation of the sample in response to a load, the detection of the resulting sample motion using OCT, and the combined effect of both processes on the final OCE elastogram. The model combines a finite element model of mechanical motion with a linear systems model of OCT image formation, using barycentric interpolation to calculate the displacement of the scattering potentials which comprise the sample. We have validated this model by comparing simulated strain elastograms of a silicone inclusion phantom against experimental data obtained using phase-sensitive compression OCE. Using this model, we have shown that the coupling between mechanical deformation and its optical detection impacts the precision of phase-sensitive OCE. We have shown how phase decorrelation, due to both local strain and to displacement, adversely affects the displacement sensitivity and leads to an optimum value of applied load. We have shown how this change in precision in turn affects the strain SNR. In principle, the model framework we present is extendable to other forms of mechanical loading and detection by OCT, providing a means of comparing and contrasting different OCE methods.