Method of Tissue Differentiation Based on Changes in Tissue Optical Properties Under Mechanical Stress Estimated with Optical Coherence Tomography

Abstract

This study highlights the possibility of contrasting the differences in tissues’ mechanical properties from Optical Coherence Tomography (OCT) data without maintaining phase stability during OCT data acquisition. The proposed method is based on the rate of attenuation coefficient changes under the mechanical pressure evaluation of the OCT data. It was shown, both on the calibrated synthetic and ex vivo biological samples, that the rate of attenuation coefficient changes observed corresponds to the sample’s mechanical properties and could be used to characterize the sample and distinguish it from other samples even if their optical properties before the pressure application are similar. This opens up the possibility to use an in vivo OCT-based system that can contrast mechanical properties without ensuring phase stability.

1. Introduction

In recent years, optical coherence elastography (OCE) has emerged as a promising tool for biological tissue evaluation in the framework of Optical Coherence Tomography (OCT) examination [1,2].

State-of-the-art optical coherence elastography methods can be divided into static and dynamic groups. Static methods are based on the application of uniaxial compression with the mechanical actuator, hence the alternative name for the approach being Compression OCE. In this case, the displacement between the unloaded and loaded states is proportional to the applied force, with Youngs’ modulus being the proportionality coefficient [2,3]. The applied force can be measured either by direct measurement [4] or through the estimation of the deformation of the auxiliary layer of a substance with known elastic properties [5,6]. In the latter case, it is assumed that the deformation in the object is proportional to the deformation in the auxiliary calibration layer, and their ratio allows one to estimate the Youngs’ modulus of the object from the known Youngs’ modulus of the auxiliary layer. To estimate the deformation, the phase changes under the load are measured [7,8]. Phase-sensitive OCT can be utilized for such measurements [9,10], which requires the high relative stability of the object and the OCT system during data acquisition, which complicates in vivo measurements [11].

The alternative broad class of OCE methods is dynamic methods, which can be roughly divided into wave-based OCE and vibrometric OCE. In wave-based OCE, a mechanical wave is created in the object by the locally applied brief force [12,13,14]. In the simplest case, the wave speed in the object $c_s$ is proportional to the shear modulus $\mu$ [15]:

$$\mu =\rho ·c_s^2,$$

(1)

where $\rho$ is the object’s density (∼1000 kg/m³ for biological tissues). Additional coefficients and terms can be applied to Equation (1) depending on the particular method of the wave excitation [16]. Youngs’ modulus $E$ can be calculated using the shear modulus $\mu$ as follows [2]:

$$E=2\mu ·(1+\nu ),$$

(2)

where $\nu$ is Poisson’s ratio, which for soft tissues at small deformations can be approximately set to 0.5. Hence, Equation (2) can be written as follows:

$$E\cong 3\mu .$$

(3)

For wave speed measurements, the Doppler shift is calculated using the phase difference between the subsequent OCT frames [13].

Vibrometric OCE is based on the measurements of local oscillations created by the transient forces applied to the tissue [17,18,19]. In this case, the investigated tissue can be described as a spring–mass damper system, in which displacement is defined by the following equation [19]:

$${\displaystyle \frac{d^{2}y(t)}{d{t}^2}}+4\pi ·f·{\displaystyle \frac{dy(t)}{dt}}+{\left(2\pi f\right)}^2·y(t)=0,$$

(4)

where $y(t)$ is displacement, $\zeta ={\displaystyle \frac{c}{2·\sqrt{m·k}}}$ is the damping ratio, $f={\displaystyle \frac{\sqrt{\frac{k}{m}}}{2\pi }}$ is the damped natural frequency, $c$ is the viscosity coefficient, $k$ is the spring constant, and $m$ is the equivalent mass. The damped natural frequency $f$ allows us to estimate the tissue’s Youngs’ modulus as follows [20]:

$$E=k·{\displaystyle \frac{L}{S}} ,$$

(5)

where $L$ and $S$ are the thickness and the contact area of the sample. The ability to accurately measure the damped natural frequency is also dependent on the phase stability of the OCT system [18].

The phase stability is essential for the majority of the state-of-the-art OCE methods, which requires additional efforts, especially if the system is intended to be used for in vivo measurements [11].

One of the advantages of the OCT setup for biological tissue evaluation is the potential to use the OCT data for extracting other tissue characteristics. One of the commonly used characteristics is the attenuation coefficient. For instance, it was widely used for the evaluation of brain tissues, including brain cancer margin assessment [21,22]. It was also used for the evaluation of burn scar tissue [23] and lymphatic nodes. In a recent study [24], it was demonstrated that the use of the elastography information with the optical attenuation coefficient extracted from the same dataset can be beneficial for diagnosing breast cancer resection margins accurately. The attenuation properties themselves may change under mechanical pressure [25,26]. Previously, it was suggested that the change in the optical properties of biological tissues may be associated with the compaction of scatters [27,28] or a change in the refractive index due to liquid extrusion [29,30]. The qualitative change in the scattering properties of the tissue was used as a diagnostic tool in [31], where the authors proposed compress the investigated tissue to better distinguish between the tumor and edema. In the present study, we propose to measure such changes quantitatively with the depth-resolved attenuation coefficient estimation method [32] and use the dependence of such changes on the applied stress to evaluate the sample. Unlike most of the state-of-the-art methods for OCT-based elastography, the proposed method does not require phase stability while collecting the OCT data and thus could potentially be applied to in vivo measurements without overcomplicating the OCT setup.

2. Materials and Methods

2.1. Samples

Several samples were evaluated with the proposed method. In the first part of the study, two-component molding silicones (Epoxy Master, Moscow, Russia) were used as a test object. After mixing both components comprising the compound, two samples with a Shore Hardness Index of 10 A and 30 A were molded using a prefabricated form.

Four samples were evaluated in the second part of the study, soft and hard plastic (polyurethane) phantoms (INO’s Biomimic TM, Québec, QC, Canada) as well as chicken skin and muscle tissue, to demonstrate the feasibility of the proposed approach for the biological tissues.

2.2. Evaluation of Samples Hardness

An REX 1600 Shore-type OO durometer (REX Gauge company, Buffalo Grove, IL, USA) was used for an independent assessment of the studied samples hardness. The cylindrical samples with a diameter of 23 mm and a thickness of 10 mm were made to measure the hardness of 10 A and 30 A silicones. For the soft plastic phantom, a cylindrical sample with a diameter of 15 mm and a thickness of 10 mm was measured. For chicken muscle tissue, measurements were carried out using a 10 mm thick rectangular sample with a lateral size of 20 × 20 mm². The axial size of the test samples was significantly larger than the length of the indenter of the measuring device (2.5 mm), which ensured the correct measurement of the sample’s hardness. Five hardness measurements were taken at different points for each test sample. The average values were then calculated.

Measurements were not performed for two samples: the hard plastic phantom and the chicken skin tissue. The hard plastic phantom is too hard to measure on the OO Shore scale as it is a hard polyurethane with a Shore Hardness Index of 84 D. The second sample, the chicken skin tissue, is too thin to measure the hardness of accurately.

2.3. OCT Setup

The study was performed with a spectral-domain OCT (SD-OCT) device developed at the Institute of Applied Physics of the Russian Academy of Sciences (Nizhny Novgorod, Russia). The device has a 20,000 A-scan/s scanning rate and performs 2D lateral scanning with a range of 1.5 × 1.5 mm² to obtain a 3D distribution of backscattered light in the polarization parallel and orthogonal to the polarization of the probing beam. The device operates at a 1 µm central wavelength and has a lateral resolution of 15 µm.

2.4. Experiment for Attenuation Coefficient Change Estimation Under Uniaxial Force

The experimental setup is shown in Figure 1. The OCT probe was fixed perpendicularly to the sample. The sample was placed on the uniaxial translation table. The axial movement of the table was used to create and adjust stress on the sample. The experiment consisted of two parts. The first part included model measurements on the silicones with the known Shore Hardness Scale. The second part was designed to emulate the compression OCT. The difference between the experiments was in the method of measuring the pressure and recording the OCT data.

In the first part, the PV-30 scales (MASSA-K, St. Petersburg, Russia) were used to measure the stress on the sample. In the second part, the layer of the auxiliary low-scattering silicone was placed between the OCT probe and the sample for calibration and stress control.

In the first experiment, the OCT image was recorded at several stress values of. In the second experiment, the OCT dataset was recorded with increasing stress on the sample, i.e., one recorded volume corresponded to the range of the applied stresses.

2.5. Auxiliary Silicon Layer Calibration

To measure the dependence of the auxiliary silicon layer thickness on the applied stress, the calibration experiment was performed. The silicon sample on the thick scattering piece of plastic was placed on the PV-30 scales (MASSA-K, St. Petersburg, Russia). The additional plastic layer was used to avoid back reflection from the metallic scales’ surface. The OCT probe on the axial translator was placed above the silicon sample, and the stress applied to the sample by the probe’s axial movement was measured by the scales. The auxiliary silicone layer thickness was measured from the OCT image of this layer. The obtained dependence of the layer thickness on the stress was approximated with the quadratic polynomial function. This function was used in the experiments to estimate the applied stress from the auxiliary silicone thickness.

2.6. Attenuation Coefficient Estimation

To quantitatively estimate the attenuation coefficient, the depth-resolved method proposed in [33] was used. This method is the extension of the method proposed in [34] and corrected for the additive noise. The attenuation coefficient was estimated from the measured OCT signal as follows:

$$\begin{gathered} \mu [i]={\displaystyle \frac{H\left[i\right]·{SNR}^{\mu }\left[i\right]}{{\left|H\left[i\right]\right|}^2·{SNR}^{\mu }\left[i\right]+1}}·{\mu }_{est}[i], \\ H\left[i\right]=1-{\displaystyle \frac{{\sum }_{j=i+1}^{\infty }N\left[j\right]}{{\sum }_{j=i+1}^{\infty }\left(I\right[j]+N[j\left]\right)}}=1-{\displaystyle \frac{\langle N\rangle ·(i_{max}-i)}{{\sum }_{j=i+1}^{i_{max}}\left(I\right[j]+N[j\left]\right)}}, \\ {SNR}^{\mu }\left[i\right]={\sum }_{x\_i,z\_i\in W}{\displaystyle \frac{{\left|{\mu }_{est}\left[i\right]\right|}^2-{\left|N^{\mu }\left[i\right]\right|}^2}{{\left|N^{\mu }\left[i\right]\right|}^2}}, \\ {N}^{\mu }\left[i\right]={\displaystyle \frac{N\left[i\right]}{2·∆·{\sum }_{j=i+1}^{\mathrm{\infty }}(I\left[j\right]+N\left[j\right])}}={\displaystyle \frac{\langle N\rangle }{2·∆·{\sum }_{j=i+1}^{i_{max}}(I\left[j\right]+N\left[j\right])}}, \\ {\mu }_{est}\left[i\right]={\displaystyle \frac{I\left[i\right]+N\left[i\right]}{2·∆·{\sum }_{j=i+1}^{\infty }\left(I\left[j\right]+N\left[j\right]\right)}}, \end{gathered}$$

(6)

where $\langle N\rangle$ is the mean noise amplitude, ${SNR}^{\mu }\left[i\right]$ is the local signal-to-noise ratio, $(I\left[j\right]+N\left[j\right])$ is the measured OCT signal, $∆$ is the pixel axial size, $i$ is an axial measurement number, $i_{max}$ is the total number of pixels in the axial direction, and ${\mu }_{est}\left[i\right]$ is an attenuation coefficient value estimated according to [34]. The estimated attenuation coefficient was averaged across the sample to decrease the noise influence on the results.

2.7. Signal Processing in Silicone Experiments

In the experiments with the molded silicones, the applied stress was estimated using the values measured by the scales, and the resulting deformation was measured as the changes in the silicone sample thickness. The attenuation coefficient was estimated according to Equation (6). To plot the dependence of the estimated attenuation coefficient on both applied stress and deformation, the attenuation coefficient and deformation values were both averaged across the B-scan. The dependencies of the measured stress applied to the object on the object’s deformations were also recorded and linearly approximated to estimate Youngs’ modulus of the $E$ silicones.

In the second set of experiments, objects with a thickness far exceeding the OCT imaging depth were evaluated, which made the direct estimation of their deformation challenging. The stress applied to the sample was estimated using the deformation of the calibration layer. The estimation was made for each A-scan. The attenuation coefficient was estimated with Equation (6). In the lateral dimension, the attenuation coefficient was averaged over 64 regions, each being 4 pixels wide. In the axial dimension, the averaging was performed in the range free from the attenuation coefficient estimation artifacts. Near the surface of the sample, such artifacts were caused by the Fresnel reflection from the sample border. Near the bottom of the visualized samples, such artifacts were introduced either by the estimated attenuation value suppression in the low SNR according to Equation (6) in the case of high-scattering biological samples or the remaining OCT signal, which was not fully attenuated at the maximal B-scan depth in the low-scattering plastic phantoms. Thus, the averaging was performed with an indent of 30 μm from the surface of the sample as well as from the maximal depth of the visualized cross-sectional distribution of the attenuation coefficient. Thus, for each B-scan, 64 pairs of attenuation coefficient-applied stress were obtained, and the scatterplot for all 512 B-scans was plotted.

2.8. Statistical Evaluations

To estimate the dependence of the attenuation coefficient on the applied stress, the linear regression was used:

$$y_{\mu }=-\alpha ·P+{\mu }_0,$$

(7)

where $y_{\mu }$ is an attenuation coefficient value, ${\mu }_0$ is the bias of the linear regression model representing the attenuation coefficient value in the absence of the stress, $\alpha$ is the linear coefficient representing the rate of the attenuation coefficient changes under the applied stress, and $P$ is the applied stress.

The confidence intervals for the regression coefficient can be calculated as follows:

$$∆\alpha =t_{a/2,n-2}·\sqrt{{\displaystyle \frac{\frac{1}{n}·\sum {(\mu -y_{\mu })}^2}{\sum (P-\overline{P})}}},$$

(8)

where $\mathsf{\Delta }\alpha$ is a confidence interval for the coefficient $\alpha$, $t_{a/2,n-2}$ is a critical value that depends on $a$ and $n$, $n$ is the number of the experimental data points, and $a$ is a parameter that determines the width of a confidence interval. In this case, $(1-a)·100\%=95\%$. The summation is performed over all data points.

The dependence of the attenuation coefficient on the deformation of the silicone layer was estimated in a similar manner:

$$y_{\mu }=-\beta ·S+{\mu }_0,$$

(9)

where $\beta$ is the coefficient describing the dependence of the rate of the attenuation coefficient changes on the deformation and $S$ is the deformation.

2.9. Parameters Mapping

The results of the experiments were plotted on a 2D scatterplot. The vertical axis shows the averaged value for the attenuation coefficient of a sample with the interval depicting a standard deviation for the measurements. The horizontal axis shows the regression coefficient calculated according to Equation (7) with the intervals calculated according to Equation (8).

3. Results

3.1. Hardness of Samples

The measurements showed that silicone 30 A has a Shore Hardness Index of 70.4 ± 2 OO, that is silicone 10 A is 47.6 ± 2 OO, that of the soft plastic phantom is 64.5 ± 2 OO, and that of the chicken muscle tissue is 30 ± 2 OO.

The hardness of the hard plastic phantom is outside the Shore OO scale (84 D). The chicken skin tissue is too thin to measure its hardness correctly.

3.2. Experiments with the Molding Silicones

In Figure 2, the dependencies of the applied stress on the sample’s deformation, the attenuation coefficient on the deformation and the applied stress for two molding silicones with Shore Hardness values of 10 A and 30 A are plotted. The linear approximation of the stress dependence on the deformation results in a Young’s modulus equal to 198 and 639 kPa for the molded silicones. The graphs presented in Figure 2 show that with the low values of applied stress <10 kPa, the dependence is nonlinear, which may correspond to the process of the silicone distribution on the imperfections of the surface. The rate of the attenuation coefficient dependence on the applied stress $\alpha$ is higher for the soft silicone (2.1 × 10⁻³ 1/(mm × kPa)) than for the hard one (1.7 × 10⁻³ 1/(mm × kPa)), while that of the rate of the attenuation coefficient changes on the deformation shows the opposite behavior, as the value for the silicone 30 A is equal to 0.99 1/mm and that for the silicone 10 A is equal to 0.44 1/mm.

3.3. Auxiliary Silicon Layer Calibration Curve

The calibration experiment results and the corresponding quadratic approximation are shown in Figure 3.

3.4. Plastic Phantomes

In Figure 4a,b, the B-scans of the two studied plastic phantoms are shown, and the correspondent cross-sectional distributions of the attenuation coefficient estimated with Equation (6) are shown in Figure 4c,d. The estimated attenuation coefficients show the linear dependencies on the stress (see Figure 4e,f), and the parameters of such fits for two evaluated phantom samples show statistically significant differences.

3.5. Biological Samples

The OCT B-scans of the chicken muscle and skin tissues are shown in Figure 5a,b, respectively. The corresponding attenuation coefficient cross-sectional distributions are shown in Figure 5c,d, respectively. The skin tissue sample shows nonlinearity at high stress. The dependence for the muscle tissue is linear. The tilts of these dependencies show statistically significant differences (see Figure 5e,f).

3.6. Map of Parameters

Figure 6 shows the map of the parameters: the average attenuation coefficient and the linear estimated rate of the attenuation coefficient changes. The confidence intervals of the six studied samples do not intersect with each other on the map of two parameters. However, the samples in pairs have a similar attenuation coefficient and cannot be distinguished by this parameter alone. The introduction of an estimate of the rate of the attenuation coefficient changes has the potential to improve the differentiation of tissues. Moreover, several samples have a similar slope parameter, which may indicate that they have similar mechanical properties.

4. Discussion

The present study demonstrates the possibility of utilizing the rate of attenuation coefficient changes under pressure as an additional characteristic for tissue differentiation. The proposed approach allows the characterization of a sample with both the attenuation coefficient and its dynamics under the applied stress from the same OCT data. All the studied samples showed a decrease in the estimated attenuation coefficient under stress, which a change in the statistical distribution of scatterers can explain. This may be a process similar to the reduction in the average refractive index difference due to the extrusion of liquids. The effect of reducing the attenuation coefficient was also observed in studies by other authors [30,35]. A linear regression model was used to estimate the change rate, but the attenuation coefficient’s actual dependence on stress may have been nonlinear. We assume that the effect of this potential nonlinearity on the rate of the change estimate would not have been significant in this work. The obtained results suggest that the rate of the attenuation coefficient changes is related to the mechanical properties of the materials: in the small data sample studied, the higher Shore Hardness corresponds to the lower rate of the attenuation coefficient changes; however, more experiments are needed to confirm this dependence. It is impossible to correctly measure Shore Hardness for chicken skin tissue because this object is too thin. The rate of change in the attenuation coefficient for chicken skin tissue is significantly higher in comparison with that for other objects. This may be due to the fact that this tissue was initially very thin, and the relative changes, with applied stress both for the thickness of the sample and for its optical properties, were much higher for this sample.

At the same time, the values of the attenuation coefficient change rates do not provide quantitative insights into the values of the mechanical properties, at least not at the given Shore Hardness. The slope parameter is 1.7 × 10⁻³ [1/(mm × kPa)] for the silicone with a Shore Hardness of 30 A and 2.1 × 10⁻³ [1/(mm × kPa)] for the silicone with a Shore Hardness of 10 A. The ratio of parameters is 1.26, while the ratio of attenuation coefficients is 1.44. The hard plastic phantom showed a value of 0.8 × 10⁻³ [1/(mm * kPa)], while the soft one showed a value of 2 × 10⁻³ [1/(mm × kPa)]. The slope parameter for soft plastic is 2.55 times greater, while the ratio of the attenuation coefficients is only 0.99. For the biological tissue, the estimated attenuation coefficient change speed rate was 2.54 × 10⁻³ [1/(mm × kPa)] for muscle tissue and 8.2 × 10⁻³ [1/(mm × kPa)] for skin tissue. The slope parameter for the skin tissue is 3.25 times greater, while the ratio of attenuation coefficients is only 0.96.

The relation of the stress and the deformation for the auxiliary silicon layer shows a nonlinear dependence, which can be explained by the absolute values of the deformations and the way the deformations were estimated. In [36], the different approaches for deformation estimations are compared. The authors concluded that the cumulative deformation estimation, carried out by the cumulative summation of the interframe deformations, leads to a better congruence with Hooke’s law. In the present study, the deformation was estimated as a difference in the current thickness of a silicone layer from its non-deformed thickness. While this method aligns with the deformation definition, it might not be optimal in real-world applications.

The nonlinear dependence was approximated with a second-order polynomial (see Figure 3). The approximation was used for stress estimation in further experiments. The exact form of the function that describes this nonlinear dependence is unknown to us. A quadratic polynomial was used for approximation because this function form provides the smallest deviation from the experimental data. The use of such an approximation cannot affect the quality of subsequent measurements.

While other elastography methods allow for the detection of displacements with subpixel accuracy, the proposed method does not measure displacements and aims to detect the differences in object’s optical property changes under the applied stress. When emulating the elastography method, the deformation of the silicone layer was assessed by changes in the thickness of the layer. This method of determining stress does not have a sub-pixel resolution, but for the purposes of this work, this is not important.

The main limitation of the proposed method compared to traditional OCT-based elastography methods is that the introduced assessment technique only allows us to distinguish between materials with different mechanical characteristics. Traditional OCE methods allow us to directly estimate the value of Young’s modulus. However, the phase stability requirements of the traditional methods are challenging to ensure in in vivo studies, while the approach proposed in this work does not require the same system stability. In addition, the new method does not require significant changes in the standard OCT system configuration.

For the effective application of the proposed method to heterogeneous media and tissues with a nonlinear characteristic of the attenuation coefficient change, a number of improvements are required. For fundamentally heterogeneous tissue, it makes sense to locally estimate the rate of the attenuation coefficient change at each specific point, which would be distributed over the image. While the utilized method for the attenuation coefficient estimations is depth-resolved, a direct implementation for the attenuation coefficient change traction might be unfeasible due to the low SNR in the resulting dependencies. Thus, some local averaging might be required, and the optimal trade-off between the spatial resolution of the method and the accuracy of the attenuation coefficient change estimation is a question for further research. Estimating the change from frame to frame will prevent problems with nonlinear changes in the attenuation coefficient since it will eliminate the linear estimate of the nonlinear function and will open up the possibility of tissue characterization with the particular type of such dependencies.

An additional problem in the study of heterogeneous tissues is the lack of understanding of how the optical properties of tissues should change under pressure, when different layers can influence each other. In [37], the changes in the attenuation coefficient in different layers of human skin was studied. In certain layers of the skin, there was a decrease in the attenuation coefficient, and in others, an increase. For this reason, the application of the method proposed in this work to heterogeneous biological tissues requires more in-depth study. A better theoretical understanding of the underlying processes affecting the optical properties is necessary to explain such phenomena.

5. Conclusions

A method for tissue behavior evaluation under the applied uniaxial stress using the OCT data were proposed. This method is based on the measurement of the attenuation coefficient change rate under the applied stress. It was shown that such dependencies largely vary for different samples depending on their mechanical properties. While the proposed method did not allow us to measure mechanical properties directly, it was shown that the measured dependencies of the attenuation coefficient are associated with such properties and largely vary for different samples. The major advantage of the proposed method is the fact that it does not require phase stability when collecting the OCT data. Thus, it can be easily implemented in in vivo studies. In addition, the proposed dependencies of the attenuation coefficient changes under the applied stress can be used along with the attenuation coefficient itself, which allows for a further differentiation of the different tissues.