1. Introduction

OCT is an interferometric optical imaging method capable of providing high-resolution cross-sectional images of tissue structures [1]. The diagnostic value of OCT, however, is frequently limited by lack of chemical specificity. PT-OCT [2] is a pump-probe extension of OCT with the potential to complement the structural images of OCT with co-registered chemical/molecular information [3,4]. To do so, the wavelength of the intensity-modulated photothermal laser (PT laser; i.e., pump) is selected at the peak absorption band of a molecule of interest (MOI; e.g., 1210 nm for lipid). In this arrangement, absorption of the PT light by MOIs generates localized modulated thermoelastic expansions and thermo-optic variations that can be sensed by demodulation of the phase of a time-lapse OCT dataset (probe). Co-registered structural information is simultaneously obtained from the amplitude of the same OCT signal. PT-OCT’s ability to produce chemically specific structural images would be particularly valuable for applications where the structural signals alone provide insufficient diagnostic information. For example, intravascular OCT of the coronary arteries can readily detect atherosclerotic lesions, but fails to reliably differentiate rupture prone plaques from more stable plaques [5]. Chemical lipid composition, assessed through spectroscopic measurement [6–8], has shown promise for refined vulnerability assessment, and may be directly investigated using PT-OCT. To date, PT-OCT has been employed to form depth-resolved maps with molecular specificity in human [9] and rabbit [10] tissues ex vivo, to detect gold-nanosphere-labelled cancer cells in vitro [11], and to visualize blood capillaries mouse ears [12] as well as melanin in zebra-fish eyes [3]. While in these studies the acquired PT-OCT signals were assessed qualitatively, results clearly showed that PT-OCT signals contain characteristic trends and attributes (e.g., amplitude, or overall shape) that directly correlate with the chemical and physical properties of the light-absorbing MOIs in tissue. Quantitative analysis of the PT-OCT signal, therefore, offers potential for obtaining depth-resolved maps of tissue composition. For instance, interesting studies on the measurement of blood oxygen saturation level in vessel phantoms [13] or of the concentration of ICG (a common dye used in ophthalmology) dissolved in water [14] have been reported. More recently, we showed possibility of estimating lipid concentration, as a key constituent of human atherosclerotic plaques, in tissue mimicking phantoms by quantitative analysis of the PT-OCT signals [15]. Quantitative PT-OCT imaging of real tissue, however, is complicated by the fact that the PT-OCT signal is influenced not only by the concentration of MOIs but also the optical, thermal, and mechanical properties of tissue. Decoupling the effects of MOI light absorption from other influence parameters requires refined understanding of the complex physics underlying the PT-OCT signal.

To date, several theoretical PT-OCT models have been proposed. In 2008, as the first theoretical model for PT-OCT, the variation of the optical path length (OPL) as a result of a change in the temperature field was modeled [2]. In this work, the solution of the heat conduction equation in presence of PT laser excitation was used to estimate the OPL change in liquid samples. Since this first model was developed in a 1-D space, it resulted in underestimation of key effects such as radial thermal diffusion and Poisson’s effect. In 2015, a more complex 3-D model for the PT-OCT signals in solid phantoms was proposed [16]. In this work, after calculating the illumination distribution, the thermal, and the mechanical stress-strain fields, the resulting OPL variation was calculated in a homogeneous sample. While this model served to improve the reconstruction of depth-resolved PT-OCT signals acquired in vivo, it assumed a homogeneous mono-layer geometry. A more refined model of PT-OCT in quasi-heterogenous multi-layered samples was presented in [17]. This model combined the individual components of the photothermal process to investigate the effects of PT laser power on the measured PT-OCT signal. A limitation of this model is that the mechanical expansion of the sample is modeled in 1-D space, along the axial direction, neglecting key elastic mechanical properties of a sample such as Poisson ratio or the mechanical stiffness of the surrounding medium. Moreover, heat flux that occurs between the layers of the sample was not considered in this model. We previously used a 1-D thermal-wave-based model to study the effects of PT laser modulation frequency on PT-OCT images [18]. Through this model, we confirmed that the amplitude of the PT-OCT signal is inversely proportional to the square root of the modulation frequency of the PT laser. We also investigated how an increase in modulation frequency of the PT laser improves the ability to detect two adjacent point-absorbers (aka. resolution). The key message of this work was that a compromise between the signal amplitude and spatial resolution should be considered to select the optimum modulation frequency. Recently, we also investigated the effect of MOI concentration on the acquired PT-OCT signal [15], revealing that in samples with heterogeneous thermal properties, the PT-OCT signals depend nonlinearly on MOI concentration. However, in this work, similar to some previous models, thermo-elastic expansion was only considered in 1-D. More recently, Veysset et al. published an interesting work on interferometric imaging of thermal expansion for temperature control in retinal laser therapy [19]. While the focus of this work is on determination of the optical and thermal parameters of multi-layered tissue via fitting of experimental data to a proposed comprehensive theoretical model, it highlights the feasibility and the need for similar comprehensive models in the field of PT-OCT. The above brief review suggests that while the proposed PT-OCT models developed by us [15,18] and others [2,16,17] offer key insights for better understanding of PT-OCT signals, they all suffer from key limiting assumptions. More rigorous theoretical models are needed to enable: 1) better understanding of the effects of system parameters and tissue opto-thermo-mechanical properties on experimental signals; 2) knowledge-based optimization of experimentation strategies; 3) guidance for developing reconstruction of depth-resolved tissue chemical composition.

In the present work, we propose a comprehensive model for prediction of the PT-OCT signal in heterogenous multi-layer samples considering the opto-thermo-mechanical properties of all in 3-D. The proposed theoretical model has a serial hierarchy with 3 blocks for predicting: the OCT and PT laser light fields, determining the induced thermal-wave field upon absorption of PT light by MOIs, and evaluating the subsequently induced thermo-mechanical expansion field due to the temperature change in the sample. For simplification, the model ignores possible coupling between the three blocks. The output of the blocks is then used to calculate the variation in the optical path length (OPL) resulting from the mechanical expansion and the temperature dependence of the refractive index and subsequently the variation of the OCT phase with time (aka, PT-OCT signal). As this model includes multiple layers in 3-D space, it allows us to survey the effects of the thermal and mechanical properties of the layers, such as Poisson ratio and Young’s modulus, on PT-OCT signals. Moreover, by considering the induced thermal field as a thermal-wave field, the model can reliably predict the effects of PT laser parameters (e.g., modulation frequency) on the PT-OCT signals. To gauge accuracy of the model, an experimental parametric study is carried out and discussed to realize the significance of the influence of PT laser power, modulation frequency, location of OCT focal plane, and thermal and mechanical boundary conditions on PT-OCT signals.

2. Theoretical PT-OCT model

The source of molecular-specific contrast in PT-OCT is the absorption of PT laser photons by MOIs. Considering the thermal properties of the sample, the absorbed modulated energy converts to modulated heat and generates thermal waves. As a result of diffusion of these thermal waves in the sample, a thermal field (aka. thermal-wave field) is formed. Consequently, the temperature of molecules inside this thermal field will modulate over time. Due to the thermal expansion the physical length of the sample will change, accompanied also by a variation of the temperature-dependent refractive index. These effects collectively cause a variation in OPL which can be sensed by a phase sensitive OCT system [20]. That is, if the intensity of the PT laser is modulated in a sinusoidal form at a specific modulation frequency, the ensuing temperature field and OPL will also modulate at the same sinusoidal frequency. Accordingly, by applying Fourier transformation (FT) to the acquired time-lapse OCT phase signal (aka. M-scan) and evaluating the resulting spectrum at the PT-laser modulation frequency, the modulation amplitude $\Delta \phi $ can be measured at each depth.

Based on the above sequence of physical phenomena, we propose a model to estimate the PT-OCT signal in samples in 3-D space using cylindrical coordinates. The flowchart of the proposed model is depicted in Fig. 1(a). The estimate model is comprised of three computing blocks in series: light field, thermal field, and stress/strain field. Each block is fed by the results of the previous block(s), system parameters, and material properties. Coupling between the blocks, while physically possible, is considered negligible and ignored. In this model, we assume that the sample is located between two fixed supports (Fig1.b). Depending on the specific sample geometry, the top support may be a constrained sample surface, a constrained glass slide on top of the sample surface, or an imaginary plane above sample surface in air. The fixed bottom support is at the bottom of sample, at a depth well beyond achievable imaging depths. Samples can be mono-layer or multi-layer. A layer refers to a homogeneous area of the sample with uniform opto-thermo-mechanical properties, such as a glass layer, or PDMS (Polydimethylsiloxane) layer. Then the volume between the fixed supports, in cylindrical coordinates, is sliced, resulting in a stack of disks. A portion of these disks will be located in the air above the sample, while the remaining disks will be within the mono- (or multi-) layer sample. Subsequently, the disks (slices) are assigned with opto-thermo-mechanical properties of the corresponding sample layers. This meshing approach ensures that the sample surface falls within the space between the two fixed supports which, in return, allows us to use the same boundary condition when solving for the stress/strain field regardless of number of layers in the sample and boundary condition on the sample surface.

 

Fig. 1. (a) Flowchart of proposed model; (b) schematic definition of heat affected zone (HAZ), thermo-mechanically affected zone (TMAZ), and the rigid zone in the sample as a function of thermal field in the sample; (c) top view of a slice and the location of elastic (HAZ + TMAZ) and rigid zones in it.

Download Full Size | PDF

In the first block, based on the optical properties of the sample and system parameters (e.g., PT laser power and modulation frequency), the PT laser irradiance in the sample is estimated in 3D. The temperature field in the sample is then determined in 3D considering the thermal properties and the light intensity distribution in the sample. The output of this part is the variation of temperature at every point in the sample over time. In the last block, the mechanical stress/strain field in the sample in response to the temperature change and as a function of the material’s mechanical properties is obtained. Eventually, the OPL variation and PT-OCT signal are calculated from the mechanical displacements and the temperature changes. The sections below depict the mathematical modeling of the three serial blocks. In this model, we solve the elasto-static strain equations, ignoring dynamic interactions between the various effects.

2.1 PT light field

To obtain a simple analytical expression, we assume Gaussian beam propagation, ignore scattering, and simply consider absorption-dominated attenuation of light, to estimate the intensity I of light in the sample as [21]:

Here, Z_f is the axial distance between a selected depth and the focal plane, L_ray is the Rayleigh range, and W_f is the beam waist at the focal plane which is calculated as:

In this equation, λ_in stands for the wavelength of the incident beam, L_f is the focal length of the objective lens, and $\; $W_col is the waist of the collimated beam incident on the objective lens. Projecting the beam into a sample with refractive index n, the Rayleigh range $\; $L_ray can be calculated as:

Through above equations, the first block of the model estimates the intensity of the PT beam in the sample. We should reiterate here that the above derivations are ignoring light scattering and assume that light attenuation within the effective imaging depth of PT-OCT is dominated by absorption of light. The rationale for this simplifying assumption is that in practice the imaging depth of PT-OCT is limited to at most the imaging depth of OCT, which is on the order of few scattering mean free paths. Within this imaging range contributions of scattering to the light field are expected to be significantly less than when approaching the transport mean free path or the diffuse regime. At these imaging depths, there would not be any OCT signal left. Moreover, in PT-OCT, the PT laser wavelength is intentionally selected at a wavelength at which the PT laser is very efficiently absorbed. These considerations motivated us to assume an absorption-dominated light field for our approximation. It should also be noted that the PT-OCT depends on the modulated (AC) portion of the thermal field and not the bulk temperature rise. This point is important because the AC portion of the temperature field is attenuated exponentially in spatial directions. That is, absorption of scattered light that is more than a thermal diffusion length (e.g., 4.3 µm at 2.5 kHz PT laser modulation) away from the OCT beam axis is not expected to contribute significantly to the measured PT-OCT signal. Light scattering will, however, create an offset error in the proposed model. A potential remedy to this limitation of our model is to convolve model predictions with an empirical spread function to account for scattering of light (as shown in Ref. [19]).

2.2 Thermal field

Absorption of photons’ energy by MOIs generates heat in the sample. The partial differential equation governing the thermal diffusion field in time-domain is [22]:

Here $\tilde{\theta }$ is the temperature in the temporal frequency domain, $\tilde{Q}$ is the heat source distribution in the frequency domain, $\tilde{I}$ is the intensity in the frequency domain and $\tilde{\sigma }$ is the complex wave number which is defined as [22,23]:

Here, $\mu = \sqrt {\frac{{2\alpha }}{\omega }} $ is the thermal diffusion length. Using the Green’s function method to solve the heat equation (Eq. (6)), the temperature field for a semi-infinite volume in frequency domain can be obtained as [22]:

In the above equation, ${\tilde{G}_0}$ is the three-dimensional Green’s function. In PT-OCT, the penetration depth of the PT light is usually shorter than the thickness of the sample, so a homogeneous sample can be considered as a semi-infinite medium. For a semi-infinite medium with homogenous thermal properties, the Green’s function is [22,24]:

In multi-layer geometries, the penetration depth of the PT light may exceed the layer thickness, and we correct for this by considering the heat exchange between layers with different thermal properties and taking the thermal effusivity at the boundary of two adjacent layers into account. This effect can be applied in the model by superposing the temperature field in absence of any thermal flux with the temperature distribution in presence of a thermal flux. Assuming perfect thermal contact, the temperature at the contact surface of two layers (e.g. layer1 and layer2) ${T_m}$ is estimated as [25]:

Here e stands for thermal effusivity. The generated thermal flux q in terms of the distance with respect to the top surface $|{{z_r} - {z_0}} |$ can then be calculated as:

Once the thermal flux at boundaries is calculated via Eq. (12), the updated temperature will be calculated with Eq. (8) [22]. Subtracting the heat flux from the thermal field produces the 3D temperature field of the multilayer geometry, $T({{\boldsymbol r},t} )$.

2.3 Mechanical stress/strain field

To estimate the thermo-elastic expansion of the sample due to temperature change, the sample is radially divided into three zones (Fig. 1(b)). The inner cylindrical zone concentric with the PT beam where most of temperature rise takes place is named the heat affected zone (HAZ). The zone outside of HAZ that is not significantly affected by the temperature change, but still mechanically contributes to thermo-elastic expansion as an impedance is named thermo-mechanically affected zone (TMAZ). These two elastic zones are surrounded by a rigid zone in which temperature rise and thermo-elastic expansions are considered negligible. Fig1.c shows the location of these three zones in a top view of a slice. To find the radii of HAZ and TMAZ, first we find the depth/slice experiencing the maximum temperature amplitude inside the sample. A 10 dB radial drop in the temperature amplitude at that depth defines the HAZ radius for all depths (Fig. 1(b)). The radius of TMAZ is defined as 10 times that of HAZ (thick-walled cylinder condition). Outside TMAZ is the rigid zone where material is not displaced. Using constitutive thermoelastic equations in the cylindrical coordinate system [26,27] and boundary conditions (B.C.) reflecting above definitions of HAZ, TMAZ, and rigid zone, the radial and axial displacements (u,w) can be found as (see Supplement 1):

The goal of above system of equation is to estimate the displacement field by meshing the sample along the depth into a stack of disks with the same radius. That is, in above equations nonuniformity in distribution of axial stress (σ_z) due to thermal gradient along the depth is ignored (B.C.4 and B.C.5). This simplifying assumption enables the model to focus on calculation of normal stresses and strains as main contributors to changes to the axial physical length of each slice (i.e., needed for estimation of OPL). In addition, the total strain along the depth is considered to be zero (B.C.6.), assuming that the distance between the top and bottom sample supports is fixed (i.e., to satisfy continuity in axial direction). The steps to solve these equations are explained in more detail in the Supplement 1.

2.4 PT-OCT signal

PT-OCT works based on measurement of the change in OPL as a function of depth and at the center of HAZ via interrogation of the time-lapse (M-mode) OCT phase. OPL is defined as the product of the refractive index and the physical length of each HAZ slice, L:

The new length of each slice L₂, is calculated from Eq. 13. Changing the temperature not only alters the length of HAZ, but also affects its refractive index:

3. Experimental methodology

3.1 PT-OCT setup

The detailed schematic of our PT-OCT system is depicted in Fig. 2. The OCT and PT light sources are a 30 mW, NIR superluminescent diode with the center wavelength of 1315 nm (+/- 75 nm at 10 dB, Exalos, Switzerland), and a diode emitting at 808 nm with 30 mW power (Thorlabs, USA), respectively. A circulator is used to direct the light between the OCT laser, the spectrometer, and the 50/50 fiber coupler. First, the output light from OCT laser is delivered to the beam splitter via the circulator. Then, the portion of the light scattered by the sample and reflected by the reference mirror and transmitted through the beam splitter towards the circulator is redirected to the spectrometer. Two galvo mirrors in the sample arm serve for raster scanning the surface of the sample with the concentric OCT and PT beams, focused by an objective lens (18 mm focal length, LSM02, Thorlabs, USA). A spectrometer comprising a 2048-pixel line scan camera (Cobra1300, Wasatch Photonics, USA) with a maximum acquisition rate of 147kHz captures the spectral interference of sample and reference arm signals. Eventually, the acquired spectra are processed on a graphics processing unit (GPU) to obtain OCT amplitude and phase. Through careful instrumentation and numerical dispersion compensation we have managed to achieve performance close to that of a shot-noise-limited system, reaching a relative displacement error of 3 nm at a SNR of 35 dB, with a system sensitivity of >100 dB. The axial and lateral resolutions were measured as 10 µm and 10 µm, respectively. The maximum optical power of OCT and PT lasers on sample surface were measured as 12 mW and 5 mW, respectively. During imaging, the PT laser intensity was modulated in form of a sinusoidal waveform from 0 mW to a set maximum power (e.g., modulated between 0 and 2 mW).

 

Fig. 2. (a) schematic of the PT-OCT setup including: superluminescent diode (SLD), optical circulator (OC), spectrometer (spec) and 2048-pixel line scan camera (LSC), photothermal laser (PT), 50:50 fiber coupler, polarization controller (PC), collimator(C), dispersion compensation block (DCB), reference mirror (RM), reflective collimator (RC), 2 degree of freedom galvo mirrors (GM), and objective lens (OL); (b) Detailed view of the sample arm.

Download Full Size | PDF

3.2 Phantoms

For evaluating the accuracy of the model, two kinds of phantoms were prepared. The first kind used PDMS (SYLGARD184, DOW, USA) as the matrix material. The two parts of PDMS (base and curing agent) were mixed in a ratio of 1:10 by weight. To enhance light scattering, 25 mg titanium dioxide (TiO₂, Sigma Aldrich, USA) was added to 1 ml of curing agent before combination with the base; the particle-curing agent mixture was then shaken in the ultrasound bath for 10 minutes to make a homogenous liquid. A dye (IR-806, Sigma Aldrich, USA) which is soluble in methanol was selected as the absorber of the PT laser light at 808 nm. Using a UV-VIS spectrometer (Shimadzu, Japan), the absorption coefficient of a solution of the dye with a concentration of 19.58 mg/ml was measured (Table.1). The dye, dissolved in methanol, was added to the PDMS mixture and mixed well by stirring. After pouring the mixture into a mold and degasification for 15 minutes, it was put on the hot plate at 70°C for 10 h to cure. Finally, the cured sample was ejected from the mold. Four phantoms with 0.22, 0.43, 1.3, and 2.2 mg/ml absorber concentrations were made. The second kind of phantom used plastisol as the matrix material. In an aluminum dish, 10 ml of liquid plastisol was heated up on a hot plate to 110°C. Then the dye solution and 25 mg of the scattering powder were added to the boiling plastisol and stirred well for 2 minutes to make a homogeneous mixture. Next, this mixture was cast into a cube mold and was cooled in air to solidify. The concentration of the dye in this sample was 5.2 mg/ml. To make a multi-layer phantom, a slice of PDMS absorber layer was sandwiched between two layers of standard glass slides. These layers were bonded together mechanically via clear tape. The thickness of each absorber layers (4 PDMS and 1 plastisol samples) was about 1 cm.

3.3 Imaging protocol

In all experiments, 1000 points over time (M-scan) at an A-line rate of 21.6kHz were captured from each point on the sample surface (duration ≃ 46 ms). The PT laser was then turned off for 460 ms before acquiring the M-scan data of the next point on the sample surface, laterally spaced from the previous point by 5 µm. The processing steps include: background subtraction, standard spectral-domain OCT data processing, forming OCT phase M-mode datasets, DC removal of M-mode phase signals by differentiating in the temporal direction, dropping initial cycles of M-mode phase datasets to only focus on the steady-state thermal response, and calculating amplitude of phase signal by applying FT transformation and locking-in to the PT laser modulation frequency (see Fig.S.1 in supplementary) [15,17,28]. The amplitude of the phase modulation at PT laser frequency can be converted to physical displacement using [12]:

3.4 Simulation and input parameters

To examine the performance of the model, developed theory was used to carry out a series of simulations under various conditions. These conditions can be divided into two categories: conditions related to PT illumination parameters, and conditions related to material properties. In the first category, the effects of PT laser power, PT laser modulation frequency, and the location of sample with respect to OCT system focal plane were studied. The second simulation category consisted of simulation of a single-layer phantom with various concentrations of absorber, a multi-layer phantom submerged in different liquids, and an absorber layer sandwiched between layers with different mechanical stiffness. These simulations were carried out under identical illumination condition. Simulated signals were processed in the same way.

The inputs to the simulations were sample properties and system parameters. Those parameters and properties that could reliably be measured (e.g., beam size) were measured and are depicted in Table 1. For parameters that could not be reliably measured or assumed, we used an optimization code to find optimum values (e.g., thermal conductivity, thermal expansion coefficient, etc.). That is, by minimizing the loss function between simulated and experimental signals from one dataset related to the single layer PDMS sample the optimized values were obtained. The criterion for loss function was the mean square error (MSE). Once the missing parameters/properties were found from the single experimental dataset, values were used for all subsequent simulations/comparisons. The optimized values are compared to literature values in Table 2. Here, the scanning range of parameters corresponds to doubling of MSE in each direction with respect to MSE of the optimal values (i.e., the global minimum).

Table 1. Input values for the simulation

View Table | View all tables in this article

Table 2. Optimized values and reported values for material

View Table | View all tables in this article

4. Results and discussion

4.1 Effect of the PT laser light field on the PT-OCT signal

To study the effects of the PT laser parameters, a sample with 2.2 mg/ml of absorbing dye concentration (μ_a= 10cm⁻¹) was defined for the simulation. To verify the simulation results, a series of experiments with matching conditions were performed on a fabricated phantom. First, to study the effect of power of the PT laser, the sample was illuminated by the PT laser at 5 different levels of power (0, 1.07, 2.47, 3.6, and 4.5 mW). In all these cases, the focal plane was aligned with the top surface of the sample, resulting in highest intensity of the PT laser on the top sample surface. Simulated PT-OCT signals and smoothed experimental signals (5-ponints running average in spatial direction) versus depth at various power levels are plotted in Fig. 3.a. Both simulated and experimental signals exhibit monotonic increase with depth. This trend is due to the cumulative nature of PT-OCT signal in a homogenous layer with the signal amplitude at each depth encoding the effect of all preceding layers. It can also be seen that PT-OCT signals exhibit a jump in amplitude at the sample surface (i.e., the Y-intercept) followed by an increase in signal values along depth. The rate of increase in signal values is higher at depths closer to surface because in these experiments and simulations the top surfaces of the sample are exposed to higher intensity of PT beam than deeper sample areas which results in larger thermo-elastic expansion of top levels. This behavior is also predicted by Eq. (1) and Eq. (8), as larger intensity of PT laser yields larger heat source Q and subsequently larger amplitude of thermal-wave fields. Given that the movement of the sample is not restricted on the top surface in this case (i.e., zero axial mechanical stress on each slice of HAZ/TMAZ), these warmer top levels tend to expand upward, so the OPL varies drastically near the top levels which contributes to the initial jumps in the PT-OCT signals at the first few depth levels of the sample (see Fig. 3(a)). Gradually, as the light penetrates deeper in the sample, the light intensity decreases due to absorption of PT light by overlayers, and the variation of temperature and OPL in deeper levels reduces. The results of Fig. 3(a) suggest good alignment of the experimental data with the developed theory. Results also demonstrate the ability of the developed theory to predict the initial signal jump in PT-OCT signals which could not be accomplished with any of the previous theoretical models.

 

Fig. 3. (a) simulated and experimental PT-OCT signals in terms of depth at various PT laser power at modulation frequency of 1000 Hz. (b) simulated and experimental PT-OCT signals versus power in various frequencies. (c) PT-OCT signals as a function of absorption coefficient (dye concentration) at various PT laser power at modulation frequency of 1000 Hz. (d) PT-OCT signals in terms of the location of focal plane related to the sample surface at modulation frequency of 500 Hz.

Download Full Size | PDF

In Fig. 3(b), the average and standard deviation of experimental PT-OCT signals within a 20μm window below the sample surface are plotted in terms of PT laser power at various modulation frequencies. As seen at each modulation frequency, increasing the power on the sample increases the amplitude of the PT-OCT signal. These experimental observations are consistent with the predicted linear increase of the model. Deviations of experimental values from linearity are likely due to noise. Such deviations are larger at lower modulation frequencies due to the pink nature of noise in PT-OCT systems (i.e., larger error bar at the lower modulation frequencies).

To understand the effect of the absorber concentration on PT-OCT signals, single layer samples with three different concentrations of absorber (0.22, 0.43, and 1.3 mg/ml) with absorption coefficients (${\mu _a}$) of 100, 200, and 600 m⁻¹ were imaged under PT illumination at various powers (0.6, 1.07, 2.98 mW). In Fig. 3(c), experimental and simulated PT-OCT signals as a function of absorption coefficient (dye concentration) are plotted. As expected, there is a proportional relation between PT-OCT signal and dye concentration. The slopes of these lines are also directly correlated to the power of the PT laser. At higher concentrations, PT light is absorbed more efficiently, leading to the generation of more heat and eventually larger temperature variations. The greater temperature leads to a greater variation in OPL and subsequently greater PT-OCT signals. Therefore, the effects of the dye concentration and the PT laser power on the PT-OCT signals are identical and cannot be directly distinguished from each other as both parameters influence the thermal energy delivered to the sample. One possible way of decoupling the effects of MOI concentration from PT laser power is to perform spectroscopic PT-OCT at dual PT wavelengths as shown before [13,40]. However, a practical downside to this approach is the dramatic decrease in imaging time which can be potentially addressed with high-speed variants of PT-OCT introduced recently to the field [41–44]. Another key point in Fig. 3(c) is the increase in size of error bars at higher concentrations which is due to the degradation of SNR caused by enhanced attenuation of PT laser at larger dye concentrations.

To study the effect of the focal plane position, its location within the sample was changed on a single layer PDMS sample with 15 mg/ml of absorber at a constant PT laser power level (4.5 mW). The amplitude of experimental and simulated PT-OCT signals at the top surface of the sample in terms of distance between top sample surface and the focal plane are plotted in Fig. 3(d). Negative values of the distance axis in this figure indicate that the sample was closer to the objective lens than the focal plane. As seen in both experimental and simulation results, the maximum signal was obtained when the top surface of the sample was located in focus at the focal plane. As the top surface moves away from the focal plane, the PT-OCT signal of the top surface drops approximately in a symmetric manner. Based on Eq. (2), at the focal plane, the beam diameter reaches a minimum, so a focused thermal field is then generated in the sample. This focused thermal field causes greater change in temperature locally at the center of the beam than a defocused beam far from the focal plane. The more change in temperature, the more variation in OPL; thus, the amplitude of the PT-OCT signal is maximum when the top surface is located at the focal plane. Regarding this inverse relation between signal amplitude and the focal plane, the intensity of the beam becomes half after displacing the sample by the Rayleigh range (150μm for our setup). As seen in Fig. 3(d), both experimental data and simulated signals confirm this drop in intensity on order of the Rayleigh range. In addition, at off-focus planes, the OCT beam is also defocused which, in return, results in deterioration of OCT phase SNR [20]. Results of Fig. 3(d) suggest that in the interpretation of PT-OCT signals, the location of focal plane relative to the sample must be considered, because a strong absorber at deeper regions can give the same signal as that of a week absorber in focus. To conclude, the optimum location for the desired depth of the sample that needs to be imaged with PT-OCT is the focal plane.

4.2 Effect of the modulation frequency on the PT-OCT signal

To survey the effect of the PT laser modulation frequency, a single layer PDMS phantom with 2.2 mg/ml concentration of absorber was used. This sample was illuminated by the PT laser at modulation frequencies of 1000, 1500, 2000, and 2500 Hz at various power levels (1.07, 2.37, and 3.6 mW). The simulated signals and average of experimental PT-OCT signals acquired in a window of 20μm just below the sample surface as a function of modulation frequency are plotted in Fig. 4(a). As seen, there is an inverse relation between the modulation frequency and the PT-OCT signal. In other words, as the modulation frequency increases, the amplitude of the thermal-wave field, and thus the PT-OCT signal, decreases. This relation can be justified by the simulated thermal fields in the sample at various frequencies. Figure 4(b) and Fig. 4(c) show the results of simulated thermal fields in the sample at the sample surface and along the depth, respectively. These plots show the temperature variation in time domain of a section passing the center of the sample. As seen, as the modulation frequency increases, the amplitude and size of the thermal field decreases, because the penetration depth of the thermal waves (i.e., thermal diffusion length) is inversely proportional to the square root of the modulation frequency (see Eq. (7)). This trend, however, does not necessarily mean that lower modulation frequencies have priority to higher ones because in practice PT-OCT systems suffer from pink noise. In other words, at lower modulation frequencies, the PT-OCT signal is larger, but the noise floor becomes larger as well. Another relevant consideration in choosing the PT laser modulation frequency is its effect on resolving adjacent MOIs. It has been shown in a previous study [18] that at higher modulation frequencies, the resolution of the PT-OCT for distinguishing two adjacent point absorber improves. This improvement in resolution occurs because at higher frequencies, the thermal diffusion length becomes shorter, preventing the interference of the generated thermal fields of the two adjacent point absorbers.

 

Fig. 4. (a) simulated and experimental PT-OCT signals in terms of modulation frequency at various PT powers. (b) simulated temporal thermal field at the top surface, and (c) along the depth at various PT laser modulation frequencies.

Download Full Size | PDF

Results of Fig. 4 show that the SNR of PT-OCT systems is a function of the PT laser modulation frequency, and a compromise needs to be considered in selecting the optimal modulation frequency of a given PT-OCT system. In previous works, the modulation frequencies were selected from few hundred [14] to few kilo hertz [2] without any optimizing or justification. Indeed, for each sample, the optimum range of modulation frequency must be obtained specifically in light of sample properties and experimental conditions.

4.3 Effect of the surrounding medium on the PT-OCT signal

The thermal and mechanical properties of the surrounding medium of an absorber affects PT-OCT signal. To experimentally observe the thermal effect of the surrounding environment, we used the plastisol sample with higher dye concentration (5.2 mg/ml concentration of absorber). Plastisol has a low thermal conductivity, which leads to a greater thermal field amplitude in the sample. To generate two different heat fluxes on the sample surface, the sample was imaged either in air or submerged in oil. The simulated and experimental signals are plotted in Fig. 5(a); both the PT-OCT signal and the initial jump are significantly larger when the sample is imaged in air. Air and oil as the surrounding media have different thermal effusivities, so the generated thermal fields in the sample are not the same in the proximity of these media. The thermal effusivity determines how much heat can transfer at the boundary of two distinct media. The thermal effusivity of the edible oil (500–700 Ws^0.5/m²K [45]) is much greater than that of air (≃6 Ws^0.5/m²K). Based on the principles of heat transfer and Eq. (11), the sample surface temperature in case of submersion in oil is lower compared to the one in air. As a result of the larger effusivity of oil, more heat is sunk from the sample, and as a result the amplitude of the generated thermal field decreases. In Eq. (8), the last part of the right-hand side stands for the thermal flux on the surface. Because the defined Green’s function in Eq. (9), as a solution for Eq. (8), has an inverse relation with the distance from the top surface which is the location of the heat flux, the temperature drop at the top surface is greater than that of deeper levels. For the cases depicted in Fig. 5(a), the amplitude is approximately reduced by a factor of 2 in presence of heat flux, and, hence, is important to consider. These observations are consistent with a finite difference validation study (see Supplement1). Considering these results, the surrounding media can penalize the PT-OCT signal. In practice, tissues consist of water and blood that can act as a powerful heat sink and deteriorate the PT signal, and the thermal effect of the surrounding media of HAZ should be considered in the analysis of PT-OCT signal. Potentially, the amount of decrease in the initial jump could be used to measure the thermal flux in the proximity of HAZ.

 

Fig. 5. simulated and experimental signals of the sample (a) in air and submerged in edible oil, (b) free sample and sandwiched sample with glass top layer. (c) simulated signals with different values of Poisson ratio.

Download Full Size | PDF

To study the effect of mechanical stiffness of the surrounding medium on the PT-OCT signal, the PDMS sample containing 2.2 mg/ml of the absorber dye was imaged once in air and once again while sandwiched between two glass slides. The simulation and the experimental PT-OCT signals in terms of depth for the free and the sandwiched samples are shown in Fig. 5(b). There is no jump (zero Y-intercept) in the signal from the sandwiched sample. Indeed, when there is no top layer, the sample tends to expand into a direction with lower reaction, so it expands upward. However, in the sample with glass top layer, the top layer with greater mechanical stiffness (Young’s modulus≃90Gpa) opposes this expansion, resulting in compression of the layer below and/or radial expansion of the material. In this scenario, the physical length in the axial direction near the junction of these layers does not change that much (unlike in the case of the free sample). Moreover, at the junction of the PDMS and the glass layers, because of heat flux, the amplitude of temperature change will be less. OPL and PT-OCT signals from the sample with glass top layer are much lower than from the free sample because maximum temperature and the physical displacement that are the two parameters in the definition of OPL are less for the sandwiched sample. In light of these results, the stiffness of the top layer can alter the signal amplitude, and when the MOI is constrained under a top layer (e.g., subsurface lipid pools in atherosclerotic plaques), this effect should be considered. Also, the change in the initial jump because of stiffness can be considered as a criterion to determine the stiffness of the top layer. For example, as a potential application, the stiffness of the cap layer in lipid-rich coronary plaques may be evaluated by analyzing the PT-OCT top surface signal jumps to get an insight into vulnerability of plaque to rupture [46–48].

Finally, to show the effect of the Poisson’s ratio, simulated signals with/without considering the Poisson effect were studied (Fig. 5(c)). As seen, the signal amplitude in the sample without considering Poisson’s ratio (ν=0) is significantly smaller than that of a sample with appropriate Poisson’s ratio (ν=0.4). Indeed, the Poisson’s effect describes how molecules of the sample tend to displace in directions perpendicular to the direction of an applied force. After heating, when the HAZ tends to expand in radial direction, because of interaction with TMAZ, a radial pressure is generated on HAZ. Consequently, because of this pressure and Poisson’s effect, material tends to expand more in the axial direction. These results suggest that neglecting Poisson’s ratio in theoretical models leads to underestimation of PT-OCT signals. The fact that we find Poisson coefficients that are consistent with literature values when optimizing the simulation parameters offers further evidence for the validity of the model.

5. Conclusion

Photothermal optical coherence tomography (PT-OCT) is a functional extension of OCT, capable of producing co-registered maps of light absorption signatures within the high-resolution tomograms of OCT. Qualitative molecular-contrast imaging with PT-OCT has been demonstrated in tissue in vivo, ex vivo, and in vitro. However, gaining quantitative insight into chemical composition of biological tissues with PT-OCT requires refined understanding of sample and system influence parameters affecting the PT-OCT signals. In this work, we propose a comprehensive and 3D theoretical model for predicting PT-OCT responses in multilayer geometries. The proposed model considers optical, thermal, and mechanical properties of samples as well as PT-OCT system parameters such as PT laser modulation parameters and location of system focal plane with respect to sample. Simulation and experimental parametric studies presented in this work demonstrate the ability of the developed model to predict the general behavior of PT-OCT signals in multi-layer phantoms with different optical, thermal, and mechanical properties. Our results also highlight how the properties of the material in the vicinity of light absorbing molecules contribute critically to the acquired PT-OCT signals. We should, however, caution the reader that the developed model is an approximation model based on simplifying assumptions such as neglecting scattering of light in the light field block, or based on a thermal model developed for solid materials rather than tissue. These assumptions together with other random and systematic sources of noise and error resulted in suboptimal alignment of presented simulation and experimental data (e.g., in Figs. 3(b) and 4(a). Nevertheless, we anticipate the proposed model to open the door for better understanding of the effects of system parameters and tissue opto-thermo-mechanical properties on experimental signals, enabling informed optimization of experimentation strategies. The model can potentially also be used as a tool for identifying parameters that are most significant in specific experiments to subsequently guide signal processing solutions for depth-resolved prediction of tissue molecular composition information.