1. 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–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–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–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 phase-sensitive 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.

2. 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.,

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

2.2. 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-noise-limited detection, and OCT SNR >> 1, s_Δϕ ≈(SNR_OCT)^−1/2 [3]. The displacement sensitivity due to optical noise, s_dO, is then

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]

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

3.1. 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–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 shot-noise-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µ_tz), 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, $\sqrt{R(z)}$, i.e. $A(z)\propto \sqrt{R(z)}$, where $\sqrt{R(z)}=\sqrt{\rho }\exp (-{\mu }_{t}z)$. The proportional change in reflected field amplitude, $\sqrt{R_p}$, from z₁ = z to z₂ = z + δz is then

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² = x² + y²) and axial coordinate z,

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

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, σ²) [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.

 

Fig. 1 Illustration of the model of image formation in OCT, and comparison with experimental results. (a) Map of the locations of scattering potentials for the simulated OCT B-scan. (b) RM(x,z), the attenuation-scaled map of the reflected field amplitude. (c) Zoomed region of the highlighted area in (a). (d) and (e) Magnitude and phase, respectively, of the complex PSF. (f) and (g) Log-scale irradiance and phase, respectively, of the simulated OCT scan formed by convolving (c) with (d)–(e). (h) Irradiance of the added shot-noise. (i) Simulated OCT B-scan obtained by multiplying (a) with (b), convolving the result with (d)–(e), and adding (h). (j) and (k) Bright-field microscopy and experimental OCT B-scan, respectively, of the physical phantom; (j) was taken before adding scatterers to the bulk space surrounding the letters, (k) was taken after. (a)–(d) are in arbitrary units. (f), (h), (i) and (k) are scaled to the SNR of the OCT irradiance.

Download Full Size | PDF

3.2. 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:

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.

 

Fig. 2 FEM simulation of a sample containing a stiff inclusion under quasi-static compression. Black lines indicate the FEM mesh (x,z), arrows indicate the magnitude and direction of the computed local displacements under uniform compressive loading from above.

Download Full Size | PDF

3.3. 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₁, r₂, r₃) of the triangle which surrounds the scattering potential, r_s = λ₁r₁ + λ₂r₂ + λ₃r₃. The real-valued coefficients λ₁, λ₂, λ₃ are called the barycentric coordinates of r_s with respect to (r₁, r₂, r₃), and fulfill the constraint that λ₁ + λ₂ + λ₃ = 1 [48],

 

Fig. 3 Computing the new location of the scattering potentials under an applied load. (a) FEM provides a mesh of the sample geometry, and the local displacements evaluated at the vertices of the mesh. (b) Delaunay triangulation re-meshes the FEM vertices using triangular elements. (c) Barycentric interpolation uses the relative locations of the scattering potentials with respect to the triangulation to obtain the locations of the scattering potentials in the loaded sample.

Download Full Size | PDF

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.

 

Fig. 4 Flowchart of the multiphysics simulation of OCE. Blue boxes denote inputs to the simulation, green boxes denote particular processes, detailed in Section 3, and red boxes denote simulation outputs.

Download Full Size | PDF

Table 1. OCE simulation inputs and parameters

View Table | View all tables in this article

4. Experimental procedure

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

4.2. 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 λ₀ = 835 nm, and a 3-dB spectral bandwidth of Δλ = 50 nm. The objective lens in the sample arm has a focal length of f₀ = 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 free-space 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 σ², 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₂ 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⁻¹ in the inclusion and ~1.18 mm⁻¹ 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⁻¹. 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-least-squares linear regression over a fitting range of 100 μm [3].

 

Fig. 5 (a) Schematic of the sample arm of our phase-sensitive compression OCE system. Schematic of representative (b) relative displacement between the reference reflector and the sample, and (c) strain A-scans, in a homogeneous region of a phantom (blue), and through a stiff inclusion embedded in a softer surrounds (red).

Download Full Size | PDF

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

5. Results

5.1. Inclusion phantom

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.

 

Fig. 6 Experimental scans of a silicone inclusion phantom compared to results of the multiphysics simulation of phase-sensitive compression OCE. (a) Experimental and (b) simulated OCT SNR images. (c) Experimental and (d) simulated relative axial displacement within the loaded phantom. (e) Experimental and (f) simulated strain elastograms in units of milli-strain (mε). (g), (h), and (i) A-scan plots of irradiance, relative displacement and strain, respectively, of the experimental scans (blue) compared to the simulation (red). A-scans are averaged over a 30 µm lateral region indicated by the blue and red dotted boxes in (a)–(f).

Download Full Size | PDF

 

Fig. 7 Regions used for comparing experiment to simulation, shown on (a) the simulated OCT image, and (b) the simulated strain elastogram. (c)–(g) Zoomed views of the 325 µm × 50 µm (x × z) regions 1–5, respectively, from the experimental and simulated OCT B-scans, and the experimental and simulated OCE strain elastograms.

Download Full Size | PDF

Table 2. Numerical comparison of experimentally acquired vs. simulated elastograms for the regions marked in Fig. 7

View Table | View all tables in this article

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

 

Fig. 8 Displacement sensitivity (s_d) (a) vs. local strain at various depths in the sample, and (b) vs. depth at selected values of strain in the sample. Blue lines are simulation results without optical noise. Red lines are simulation results with optical noise and attenuation. Black lines are displacement sensitivity values calculated at zero strain, assuming only optical noise.

Download Full Size | PDF

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.

5.3. Strain signal-to-noise ratio

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.

 

Fig. 9 Strain SNR (SNR_ε) (a) vs. local strain at various depths in the sample, and (b) vs. depth at selected values of strain in the sample. Blue lines are simulation results without optical noise. Red lines are simulation results including optical noise and attenuation.

Download Full Size | PDF

6. 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. Zaitsev et al. [28, 29] described a method of simulating OCE images by generating a map of random values, filtered such that the correlation resembled an actual OCT image. Sample deformation was simulated using a simple analytical expression of axial and lateral displacements given an applied strain. Fu et al. [15] used a more rigorous optical simulation of OCT in their model of OCE, generating a 2-D speckle map using a random distribution of scattering particles and a linear systems model of OCT [19]; however, their study of mechanical deformation only considered homogenous samples. Chan et al. [26] used both a linear systems model of OCT and a FEM model of mechanical deformation; however, they simulated a time-domain OCT system, and implemented the displacement of the scattering potentials using a finely spaced grid. This implementation limits the resolution of the scatterer locations to the size of the grid spacing. Finer resolution of the scatterer locations in the optical simulation of the loaded sample requires finer grid spacing, with a corresponding increase in the computational memory requirements, reducing the practicality of this method to simulating larger, or 3-D, samples.

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$\Delta \varphi (\delta z)=2k\delta z+2{\tan }^{-1}\left[\left(2z_0\delta z\right)/\left(f^2-f\delta z\right)\right]$, approximated in Eq. (1) by $\Delta \varphi (\delta z)\approx 2k\delta z.$ It is clear that the error in this approximation will increase as δz increases, contributing to what we have termed, in this study, translation-induced decorrelation noise. In addition, when we consider the case of multiple scattering particles within the OCT resolution volume, displacement of the scatterers will also lead to a shift or translation in the speckle pattern, which implies that the phase difference is calculated between phase realizations comprising partially decorrelated areas of non-overlapping speckles. Similarly, strain will not only cause the mean location of the particles to displace, but will also alter the distance between the particles. In this case, the error caused by using the approximation in Eq. (1) will increase with the amount of strain, leading to what we have termed strain-induced decorrelation noise.

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.

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