1. Introduction

Controlling the temperature rise in tissue during laser therapy is fundamental to many medical procedures, from thermal ablation of tumors to dermatological and retinal treatments. Retinal photocoagulation, for instance, has been a standard procedure for decades to treat a number of retinal pathologies, such as proliferative diabetic retinopathy [1], macular edema [2], central serous chorioretinopathy [3], and others [4]. Extent of heating and the induced therapeutic effect depend on the temperature course (magnitude and duration), which, in turn, depends on the laser parameters (wavelength, power, duration, spot size, etc.) and the tissue properties (optical absorption and scattering, heat capacity, and conductivity). Recent studies have suggested that the therapeutic window, typically quantified by the Arrhenius integral Ω, of the non-damaging hyperthermia is narrower (0.1 < Ω < 1) than in conventional photocoagulation [5–7], and thus such treatments require higher precision in the temperature control.

In the retina, light absorption is primarily governed by the pigmentation in the retinal pigment epithelium (RPE) and choroid, associated with an optical absorption that is proportional to the melanin concentration [8]. Measuring the optical absorption in the retina (and in most tissues) in vivo is, however, particularly challenging since there is no direct optical access to the transmitted light, and it must therefore be done in reflection. Furthermore, the melanin concentration can vary locally within the retina and from patient to patient by a factor of two to four [9]. These variations, if not accounted for, can result in significant temperature differences and thus undesired lesions or sub-therapeutic exposures. To ensure the laser exposure be within the non-damaging therapeutic range (0.1 < Ω < 1), one can estimate using previously published finite-element thermal models of the retina [8,10] that the absorption coefficient must be determined with accuracy no worse than ±20% (see the Non-damaging treatment range section in the Supplemental Document and corresponding Dataset 1 [11]).

To assess the retinal absorption in reflection geometry, one can measure the tissue response to laser heating. A few approaches have been developed along these lines. One of them, an optoacoustic method, is based on the detection of the pressure waves at the cornea, generated by the expansion of melanosomes during repeated nanosecond laser exposures [12–15]. During the thermal and acoustic confinement, the detected pressure is proportional to the absorbed energy density, related to the absorption coefficient, and the Grüneisen parameter [16]. Precision of the temperature measurement in this approach has been limited to a few degrees. Further, the signal is sensitive to the position of the piezo sensor on the cornea and hence the measurements require a careful calibration to convert the measured pressure waves into temperature. Another approach is based on detecting the changes in back-scattered light due to microbubbles formed by vaporization, which lead to the RPE cell death [17]. However, due to associated tissue damage, this technique is not suitable for non-damaging therapies. Optical coherence tomography (OCT) has also being explored as a non-contact, high spatial (3D) resolution imaging method for in-situ monitoring of the tissue deformation [10,18–21]. In particular, phase-sensitive OCT [22] allows tracking of the nanometer-scale optical path length changes, and has been proposed for precise temperature control during photothermal therapy [18]. Yet, to date, experimental investigations have not directly and quantitatively related the phase signal to thermo-mechanical changes in tissues, from which the temperature could be derived. Instead, temperature profiles have been separately estimated via the optoacoustic method mentioned above [18] or analytically, using an a priori value for the absorption coefficient [19,20]. The calculated temperature profiles were then qualitatively compared with the phase profiles recorded via OCT, without any direct validation. Recent advances in pOCT demonstrated capabilities of detecting nanometer-scale deformations with sub-ms temporal resolution [23,24], which should enable rapid and precise optical thermometry.

In this study, we present a methodology to precisely determine the optical (absorption coefficient) and thermal (heat conductivity and coefficient of thermal expansion) properties of layered materials using line-field phase-resolved optical coherence tomography (pOCT). In this approach, we monitored the vertical displacements and the changes in the optical path lengths (ΔOPLs) of reflective and scattering layers with nanometer precision in a non-damaging regime. We developed an analytical model to the coupled thermo-elastic problem and solved it in the Hankel-Laplace transformed domain. This approach significantly increases the computation speed compared to finite element simulations, and therefore allows faster fitting of the material properties to match the spatial and temporal aspects of the experimental ΔOPLs.

2. Materials and methods

2.1 Optical setup

The optical layout (Fig. 1) followed the design by Pandiyan et al. [25] for high-speed line-scan OCT imaging. We here summarize the main components of the design; further details on the optical system can be found in Ref. [25]. The illumination is provided by a fiber-coupled super-luminescent diode (M-T-850-HP-I, Superlum, 9-mW power) with a center wavelength of 840 nm and a bandwidth of 130 nm, providing an axial resolution of 2.4 μm in air. After collimation with a reflective collimator (RC) to a diameter of 1.7 mm (FWHM of the Gaussian intensity profile), the beam is shaped using a cylindrical lens (CL1) to yield a line-illumination at the tissue phantom, which is held around its periphery and whose front and back surfaces are in contact with air.

 

Fig. 1. (a) Top-view of the optical setup including the line-field spectral-domain OCT (red path) and the laser heating (blue path). SLD: superluminescent diode, RC: reflective collimator, CL achromatic cylindrical lens doublet, L: achromatic lens doublet, BS: beam splitter (70T:30R), M: mirror, GS: galvo-scanning mirror, DFM: deformable mirror, DCM: (long-pass) dichroic mirror, TP: tissue phantom, LD: laser diode, BPF: band pass filter, P: pupil plane, I: image plane. The inset shows the overlap between the round heating beam spot and the OCT line illumination. The focal lengths of the lenses are: CL1 = CL1R = CL3 = 75 mm, L1 = L1R = L8 = 250 mm, L2 = L2R = 150 mm, L3 = L3R = 125 mm, L4 = L4R = CL2 = L2’ = 200 mm, L5 = L5R = 350 mm, L6 = L6R = L9 = 100 mm, L7 = 10 mm, L1’ = L3’ = 300 mm. The insets show the overlap between the heating beam and the OCT illumination field, and the tissue phantom assembly. (b) Simplified top view of the optical setup showing the illumination beam (solid beam) and two imaging rays in the detection path (red line). (c) Corresponding side view illustrating, with B, the anamorphic configuration based on cylindrical lenses (CL1–3).

Download Full Size | PDF

After CL1, the beam is split between the reference and sample arms with a 30:70 (R:T) beam splitter. In the sample arm, the beam, after going through magnification telescopes, is focused along the y direction to yield a line field of 0.8-mm in length (FWHM of the Gaussian intensity profile) in the x direction for an approximated focal spot size of 10 μm (beam waist diameter) and a Rayleigh range of 95 μm. A galvo-scanner (8310K Series Galvanometer, Cambridge Technology) is used to steer the beam laterally (y direction) at the sample surface and a deformable mirror (DM69-15, ALPAO) is used for fine focus adjustment and for the future adaptive optics implementation. The entrance pupil plane (P), the galvo-scanner, and the deformable mirror are optically conjugated using afocal telescopes (L1 + L2, L3 + L4).

The detection path, shown in two views in Fig. 1(b),(c) follows an anamorphic imaging configuration by using two cylindrical lenses to allow independent control of magnification in the x and y dimensions, and thus independent control of the spatial and spectral resolutions. The OCT spectrometer is constructed using a 600 line/mm grating (Wasatch Photonics) and a high-speed camera (Phantom v641, 10-μm pixel size) for a resolution of 512${\times} $ 768 pixels (spatial ${\times} $ spectral). B-scans, with a static illumination field, are recorded for 500 ms with a frame rate of 10 kHz and 98-µs exposure times for a full imaging field of view of 2.3 mm laterally and 1.0 mm in depth, one camera pixel corresponding to 4.5 µm laterally and 2.7 µm in depth. Flat glass windows are also added to the reference arm (not shown in schematic) for coarse chromatic dispersion compensation between the reference and sample arms.

For laser heating, a blue diode beam (450-nm wavelength, 1-ms duration, 260-mW peak power) is coupled into a 200-µm multimode fiber. The fiber output is imaged onto the absorbing layer with a 400-µm diameter spot of uniform intensity. An iris is introduced at an intermediate image plane to sharpen the beam image edge and adjust the spot size. The beam is coupled into the OCT path using a long pass dichroic mirror and manually steered so that the heating beam diameter overlaps with the line OCT illumination. The pulse duration and the delay for the OCT camera acquisition are controlled by a pulse/delay generator.

2.2 Tissue phantom fabrication and materials properties

A five-layer tissue phantom, emulating a simplified retina, was fabricated by spin coating several layers of poly(dimethysiloxane) (PDMS) to different thicknesses on top of a glass substrate. PDMS is a common tissue simulant material and has been used for retina phantom assemblies [26]. PDMS layers were fabricated using two-part, 10:1 mixing ratio, Sylgard 184 kits from Dow Corning, on top of a 1-inch borosilicate glass substrate (N°2 coverslip, ChemGlass) (Layer 5). The spin-coating recipes are given in Table 1. Layers were assembled in reverse layer label order, Layer 1, the top surface, being the last layer coated (see Fig. 1(a) inset). After spin-coating, samples were degassed in a vacuum desiccator for 30 minutes. All layers were then cured for 48 hours at room temperature before the next coating step. Curing at room temperature was chosen so not to degrade the dye, which visibly deteriorated at temperatures above 40°C. Since the elastic modulus is dependent on the curing temperature [27], the same curing temperature was applied to all layers to simplify the mechanical model. To make the top, scattering PDMS layer (Layer 1), titanium particles (2.5wt%, 21-nm nominal diameter, Sigma Aldrich) were added to the base PDMS, previously diluted in hexane (3:1 weight ratio 3:1) to reduce the viscosity and hence facilitate the nano-powder mixing. The hexane was left to fully evaporate during stirring and regular sonication of the mixture for 24 hours prior to mixing with the curing agent and spin-coating. Layer 3, which aims to simulate the thin scattering and absorbing RPE layer, was fabricated using a PDMS/dye/titanium oxide mixture. First, a saturated solution of hexane with Oil Red O dye, strong absorber at 450 nm, was prepared. The solution was then mixed with the PDMS base (1:1 weight ratio) to obtain a low viscosity mixture. Lastly, the titanium oxide nano-powder was added (2.5 wt% with respect to PDMS base) to the mixture, which was left under stirring and regular sonication for 24 hours. The mixture was then spin-coated at high speed.

Table 1. Tissue phantom layers and spin-coating recipes.

View Table | View all tables in this article

To measure the absorption coefficient of the dye-doped PDMS layer, a similar mixture but without the titanium particles (to avoid scattering), was prepared and spin-coated onto a glass substrate following the same recipe as for Layer 3. Light absorption in this layer, measured using the same 450-nm diode laser as in the heating experiments, was 8%. To calculate the absorption coefficient μ_a, the thickness of that layer was measured by a profilometer. For simplicity the mechanical and thermal properties of all PDMS layers were considered identical. All material properties used in the modeling and references are listed in Table 2.

Table 2. Non-fitted material parameters used in model. ρ: density, E: Young’s modulus, α_TE: coefficient of thermal expansion, ν: Poisson’s ratio, κ: coefficient of thermal conductivity, c_p: specific heat capacity, n_RT: refractive index at room temperature, α_TO: thermo-optic coefficient, μ_a: coefficient of optical absorption. References are indicated next to values.

View Table | View all tables in this article

2.3 OCT signal processing

OCT image reconstruction included the standard steps of background subtraction (of the sample and reference arms), wavenumber resampling, and time referencing to cancel the local arbitrary phase [25]. The remaining dispersion, not compensated by the glass windows, is measured by placing a mirror instead of the tissue phantom. It is then numerically compensated during the OCT reconstruction [31]. We measured a system sensitivity of 75 dB by placing a mirror at the image plane in the sample arm. The sensitivity is comparable to the system developed by Pandiyan et al. [25].

In polar coordinates (r,z), the spatially static B-scan OCT signal can be described as a complex number:

The phase difference $\mathrm{\Delta }\varphi $ between two pixel coordinates (r,z)₁ and (r,z)₂ as a function of time is calculated by taking the argument of the product of the complex and complex conjugate numbers of the OCT signals at these locations:

The OCT intensity image of the top 100 μm of the phantom is shown in Fig. 2(a). The signal-to-noise ratio (SNR), which is directly related to phase sensitivity [32], is about 20–30 dB for the scattering layers of the tissue phantom. Because the top surface is reflective, while the other layers are scattering, it shows the highest SNR (maximum value ∼30 dB) and therefore would yield the lowest phase noise. The change in phase, Δφ, can be converted to change in single-pass optical path length, ΔOPL, by:

 

Fig. 2. (a) B-scan OCT image showing the sample’s top layers. The OCT signal strength decreases with depth as the light is scattered by upper layers. (b) Changes in optical path lengths ($\overline {\textrm{AA}} $, $\overline {\textrm{AB}} $, $\overline {\textrm{AC}} $) measured at the center of the heating beam (r = 0) for the three planes A, B, and C with corresponding depths z = 0, z = 25 μm, z = 60 μm. The planes are indicated by colored arrows in (a). The ΔOPLs were measured following heating by a single laser pulse at t = 0. The phase reference point (r_ref = 505 μm, z = 0) is indicated by a yellow arrow in (a). The low SNR signal that would be obtained for $\overline {\textrm{AC}} $ using an axial reference is shown in grey in C. (c) Zoom in the first 6 ms showing gradual decrease of the OPLs during laser heating (1-ms duration).

Download Full Size | PDF

2.4 Solutions to the axisymmetric thermo-mechanical problem

The transient response of the tissue phantom to laser heating can be solved using a thermo-mechanical model for an elastic medium in cylindrical coordinates [33]. The problem can be described using the governing partial differential equations of elasticity, the generalized thermo-elastic Hooke’s law including thermal stress contributions, and the heat diffusion equation [34]. Details regarding the problem formulation can be found in the Supplemental Document. We modeled the laser heating as an external heat source flux ${q_e}$ in the absorbing layer (layer 3), using the Beer-Lambert law, where the heat flux is:

The beam spatio-temporal profile was modeled as a square pulse in time and a square top-hat profile in space convolved with a Gaussian kernel to simulate blur and defocus [35]:

To solve for the displacements and temperature fields, which are necessary to calculate the changes in optical paths (ΔOPL), we used the state space approach applied to layered media, introduced by Bufler [36] and Bahar [37]. This approach consists in reducing the partial differential equations of coupled thermo-elasticity using integral transforms to a state space equation, taking advantage of the planar layered structure. The relationship between states vectors can then be described by a transfer matrix [37,38] or a stiffness matrix [39] obtained by solving the state space equation. This method has been extensively used in macro-scale multi-layer problems in geophysics [40,41], with applications to seismology and nuclear waste management, but never at the microscale for biomaterials. The thermo-elastic problem can be analytically solved, in the transformed domain, for planar layers in a semi-infinite half space, which makes it particularly appropriate for the retina, as long as we assume a heat source radius being relatively small compared to the retina radius of curvature. In comparison, finite element methods (FEM) treat semi-infinite space by significantly extending the computational domain, therefore increasing the computation time. Moreover, for an axisymmetric configuration, inversing the solution from the transformed domain requires integration over two dimensions (the radial and temporal dimensions), whereas FEM perform integrations over three dimensions (the radial, axial, and temporal dimensions), which increases the computation time. The increased computation efficiency of the state space approach has been discussed in Refs. [42,43] suggesting a >10 × acceleration. Validation of the method and efficiency considerations can be found in the Supplemental Document (Stiffness matrix method (SMM) validation section and corresponding Dataset 2 [44]). Rapid computing of the solution is critical for fitting convergence, fitting quality evaluation via bootstrapping, and precision assessment via Monte-Carlo simulations, all of which involve calculations of many solutions. We further note that modeling a volumetric heat source is possible using the stiffness matrix method [34] but at the expense of longer computation time, which is acceptable after material properties are fitted.

We followed the formulation by Ai et al. [34] to derive the state vector ${[{{{\tilde{u}}_r}^1,{{\tilde{u}}_z}^0,\; {{\tilde{\theta }}^0}} ]^T}$, where ${\tilde{u}_r}^1$ and ${\tilde{u}_z}^0$ are the 1^st and 0^th-order Hankel-Laplace (HL) transforms of the normal displacement components in the r and z directions, and ${\tilde{\theta }^0}$ is the 0^th-order HL transform of the temperature change. Numerical integration was then performed to find ${u_r}$, ${u_z}$, and $\theta $. The derivation of the state vector and the inverse numerical transform procedures are described in the Supplemental Document (Transfer and stiffness matrix derivation section).

2.5 Arrhenius integral for determination of the therapeutic window

Thermal damage in retinal cells during transient hyperthermia has been approximated following the Arrhenius model, which assumes a first-order kinetics in denaturation of the critical molecular component, the ‘weak link’. Its concentration D(t) decreases at the rate [8]:

Values for A (1.6×10⁵⁵ s^-1) and E* (340 kJ/mol) have been calibrated and characterized in previous studies so that $\mathrm{\Omega } = 1$ corresponds to the RPE damage threshold [8]. No cellular response was observed for $\mathrm{\Omega } < 0.1$ thus the non-damaging therapeutic window was defined as $0.1 < \mathrm{\Omega } < 1$ [7]. The Arrhenius integral can be directly calculated using the temperature course at any chosen point in the tissue.

3. Results

3.1 High-speed imaging of laser heating

The phase difference can be measured between any arbitrary pair of points. However, because we are interested in measuring the absolute vertical displacements of the layers, it is prudent to use a reference point whose phase remains largely unchanged following laser heating. Here, we have chosen a point at the surface of the sample (plane A, z = 0) far enough (r_ref= 505 μm) laterally to be least affected by the thermal expansion while maintaining a high SNR (Fig. 2(a)). Selection of such a reference point is possible since, with the line-field OCT, we have simultaneous access to the phase at all radial positions in a common-path configuration (in the sample arm only) [45]. In contrast, point-scan OCT only allows common-path phase referencing along the z axis (axial referencing) at a single radial position.

The complex values of the OCT signal at the three planes were averaged over 3 pixels in depth around the plane z-coordinates (± 1 pixel) before computing the phase values $\varphi ({r,{z_A},t} )$, $\varphi ({r,{z_B},t} )$, $\varphi ({r,{z_C},t} )$. The phase reference $\varphi ({{r_{ref}},{z_A},t} )$ was obtained by averaging the complex values of the OCT signal at the plane A z-coordinate and a radius offset of ${r_{ref}}$ = 505 µm (± 1 pixel in r and z). The phase at ${r_{ref}}$ = -505 µm was also calculated to retrieve and cancel the surface tilt angle in the r-z plane.

The optical path difference was tracked as a function of time for the three planes of interest: the free surface (plane A, z_A = 0), an intermediate plane in layer 1 (plane B, z_B= 25 μm), and at the absorbing layer (plane C, z_C = 60 μm). For those planes, the ΔOPLs can be expressed as:

Figure 2(b) shows the $\overline {\textrm{AA}} $, $\overline {\textrm{AB}} $, and $\overline {\textrm{AC}} $ at the center of the heating beam (r = 0). The experimental OPL noise, 5.1 nm at r = 0, z = 0, is close to the theoretical value, 4.0 nm for a phase difference calculated between points with SNR values of 30 dB and 22 dB (lateral offset point) [32].

Assuming for simplicity of interpretation that the reference phase is unchanged, $\overline {\textrm{AA}} $ is directly opposite to the top surface displacement. Upon laser heating, the temperature of layer 3 increases and heat diffuses to adjacent layers. Heating leads to thermal expansion, so the layers deform and become thicker. It causes the top layer to move up (toward negative z) and $\overline {\textrm{AA}} $ to decrease by the vertical displacement of the plane A (${u_z}({r,\; {z_A} = 0} )$). Meanwhile, as the temperature increases by θ, the index of refraction of PDMS decreases by $\Delta n$, with $\Delta n = {\alpha _{TO}}\theta $, where ${\alpha _{TO}}$, the thermo-optic coefficient of PDMS, is negative. Therefore, as plane B moves up, the optical path $\overline {\textrm{AB}} $ decreases by the additive effects of a shorter physical distance and a lower index of refraction. Since plane A also moves up and air is replaced by the polymer over the displacement thickness, the optical path $\overline {\textrm{AB}} $ decrease is slightly reduced by this effect. Same goes for $\overline {\textrm{AC}} $. Because the thermal expansion (${\alpha _{TE}})\; $ and thermo-optic (${\alpha _{TO}}$) coefficients are intimately related, their combined effects resulted in $\overline {\textrm{AA}} $, $\overline {\textrm{AB}} $, and $\overline {\textrm{AC}} $ being of similar amplitude (Fig. 2(b)). The simplified expressions of $\overline {\textrm{AA}} $, $\overline {\textrm{AB}} $, and $\overline {\textrm{AC}} $ are:

3.2 Retrieval of the material properties

We assumed that the properties of the glass substrate are well known and, given the separation between the top of the assembly and the glass, that the contribution from the glass to the response was very small. We considered for simplicity that all PDMS layers have the same isotropic thermal and mechanical properties, including the layers containing titanium oxide particles, given their low weight fraction (2.5 wt%). The density ρ, the Poisson’s ratio ν, and the specific heat capacity c_p are all well characterized for PDMS [28,29] and are considered fixed parameters. Similarly, we assumed the index of refraction n at room temperature (23°C) and the thermo-optic coefficient α_TO [30] to be known. On the other hand, the coefficient of thermal expansion α_TE strongly depends on the curing temperature, which leads to uncertainty in its value [28], and a range of coefficient of thermal conductivity κ can be found in the literature. We therefore considered α_TE and κ as fitting parameters. Although the Young’s modulus E depends on the curing temperature [27], it has little effect on the state vector (see Supplemental Document, Elastic modulus effects section). Lastly, the absorption coefficient ${\mu _a}$ was fitted through the deposited heat magnitude ${q_0}$ with ${q_0} \propto 1 - {e^{ - {\mu _a}l}}$, where l is the absorbing layer thickness. We further assumed that these properties are temperature independent, given the small temperature range of the study (<10°C). The Gaussian blur to the ideal disk (top-hat) intensity profile (Eq. (5)) [35] added two geometrical fitting parameters to the model (width of the blur σ_defocus and top-hat profile radius a_laser). This brought the total number of fitting parameters to five: α_TE, q₀, κ, σ_defocus, and a_laser.

Fitting of the parameters was performed in two steps. First, these properties were fitted by minimizing the sum of the squares of the residuals between the model and the experimental data $\overline {\textrm{AA}} $, which showed the highest SNR, using the exact ΔOPL expression (equation S2). The laser geometrical parameters, σ_defocus and a_laser, were mainly fitted through the spatial shape of the response, and the coefficient of thermal conductivity was mainly fitted through the temporal decay of the response. Importantly, however, $\overline {\textrm{AA}} $ scales linearly with both α_TE and q₀ so fitting of $\overline {\textrm{AA}} $ did not allow decorrelation between α_TE and q₀, therefore the product ${\alpha _{TE}} \times {q_0}$ was fitted through the magnitude of the response. Second, to find α_TE and q₀, we used the experimental data $\overline {\textrm{AC}} $, which is only proportional to q₀. The fitted parameters found in the first step were used to solve for $\overline {\textrm{AC}} $, and then amplitude of $\overline {\textrm{AC}} $ was fitted to determine the absolute values of α_TE and q₀, as described in the Supplemental Document (Material property fitting procedure section). The comparison between the fitted model and the data is shown in Fig. 3 for $\overline {\textrm{AA}} $ and $\overline {\textrm{AC}} $. The data are well reproduced in both the spatial (radial) and temporal dimensions, thus suggesting high fidelity of the model. The uniformity of the residuals (see Supplemental Document) over the entire data set also attests for the accuracy of the model and the absence of systematic errors. The resulting fitted values are listed in Table 3. We note a good agreement between the fitted and initial values for α_TE (286 vs 310 K^-1), κ (0.175 vs 0.17 W.m^-1.K^-1), and μ_a (20.3 mm^-1 vs 19.9 mm^-1). The uncertainties in parameter determination were obtained by fitting bootstrapped data, resampled over a total of 200 iterations, and are listed in Table 3. The obtained uncertainties are an order of magnitude lower than for the starting values of α_TE and κ, hence demonstrating the high precision of the parameter determination using this method. We also note that the obtained precision on μ_a (2.5%) is much higher than the required accuracy for therapeutic titration (20%). We finally verified the low correlation between the parameters (see Parameter correlations section in Supplemental Document) obtained by repeated fitting of the bootstrapped data.

 

Fig. 3. (a) $\overline {\textrm{AA}} $ and (b) $\overline {\textrm{AC}} $ experimental data (left halves) compared to the fitted model (right halves) in space (radial) and time. The temporal profile of the response at the center of the heating beam (r = 0) and its radial profile shortly after the heating pulse was turned off (t = 1.5 ms) are shown on the right and below the maps. The model was fitted on the rolling average of the raw experimental data.

Download Full Size | PDF

Table 3. Non-fitted and fitted PMDS material parameters. α_TE: coefficient of thermal expansion, κ: coefficient of thermal conductivity, μ_a: coefficient of optical absorption.

View Table | View all tables in this article

3.4 Temperature estimation and precision

Since we are ultimately interested in the temperature, the best-fit parameters can be used to compute the temperature increase θ at any point in space and time. For non-damaging retinal therapy, the temperature rise at the RPE, represented here by the absorbing PDMS layer, is of utmost importance. Figure 4(a) shows the temporal course of the temperature at that layer. Temperature increased during the 1 ms laser pulse, and then decreased slowly due to heat diffusion. Figure 4(b) shows two snapshots of the temperature distribution at the end of the heating pulse (t = 1 ms), when heat is relatively confined, and after 4 ms, when heat has diffused radially and axially. To estimate the temperature precision of the calculations, we ran Monte-Carlo simulations [46] using the fitted parameters with their uncertainty over 1000 iterations to compute the peak temperature (t = 1 ms) at the center of the beam at the absorbing layer (r = 0, z_C = 60 μm), and found a temperature uncertainty of ±0.25°C. These simulations at a single point in time and space took about 20 seconds. For comparison, running the same number of iterations using the COMSOL model would take approximately 6 hours.

 

Fig. 4. (a) Radial and temporal distribution of temperature at the absorbing layer (z_C = 60 μm). The temporal course at the center of the heating beam (r = 0) calculated with the best-fit parameters (black line) is shown on the right. Monte-Carlo simulations (n = 50) considering parameter uncertainties are shown in orange to illustrate the uncertainty in the temperature estimation. (b) Spatial distributions of temperature at 1 and 4 ms, showing heat diffusion into the sample.

Download Full Size | PDF

4. Discussion

Determination of the temporal and spatial distribution of the temperature rise in tissue is critical for evaluating the safe therapeutic window and the extent of potential thermal damage. Both are typically assessed using the Arrhenius equation, which describes the first order kinetics of the thermal denaturation. The therapeutic window for non-damaging hyperthermia corresponds to the Arrhenius integral value between 0.1 and 1, where values above 1 correspond to permanent damage. With the 2.5% precision of the absorption coefficient determination in our approach, the uncertainty in $\mathrm{\Omega }$ is about 0.17, which is well within the therapeutic window. Uncertainty of 0.50 in Arrhenius integral, approaching the boundary of the therapeutic window corresponds to 2.5 times lower precision in determination of the absorption coefficient (i.e., 6.25%). At the current laser settings, this would translate to a temperature rise uncertainty of about 0.65°C. Precision of this method exceeds the other methodologies for temperature-controlled retinal treatment and is more than sufficient for titration the non-damaging thermal therapy. In the current example, the absorption coefficient determination and the temperature calculations are performed way below the damage threshold (Ω < 10⁻³), which enables repeated measurements (multiple heating pulses) for further data averaging.

We are working on adapting this approach to in-vivo measurements and hope to achieve a similar level of the temperature precision despite the potential difficulties associated with tissue movements and reduced OCT signal strength, which directly affects the phase noise. Here, with the tissue phantom, we only considered the phase information from two planes (A and C), as they enable sufficient precision in determination of the material parameters. However, more data points from additional planes (e.g., plane B) are available for fitting, which would increase the precision and allow additional fitting parameters. On one hand, retina offers many scattering planes from which phase can be extracted, while on the other hand, we expect a larger number of variables for the retina. It remains to be seen whether the increased number of experimental data points provided by the additional layers will compensate for the increased number of unknown parameters.

Retinal thermal expansion coefficient is often assumed to be equal to that of water [19], but a recent study suggested otherwise and estimated it to be four times larger [18]. We can also expect this coefficient to be anisotropic (albeit transverse isotropic) for the retina, given its anisotropic cellular structure. The same might apply to the coefficient of thermal conductivity and, for the same reason, layers could exhibit different Poisson’s ratios. Conversely, the index of refraction of the retina is relatively uniform within 1% across the layers and has been reasonably determined [47]. Questions remain, however, regarding the thermo-optic coefficient (${\alpha _{TO}}$), which might more strongly depend on the solid content [48,49]. It also has been shown that the elastic modulus varies throughout the retina [21], and is two-to-three orders of magnitude lower than for PDMS (kPa instead of MPa). We therefore expect it to play a more significant role than in the tissue phantom. In addition, the present model assumes a uniform laser irradiance and a uniform absorption coefficient within the laser spot. Additional contributions from metabolic heat generation and convection due to blood perfusion can be implemented in the model, although the metabolic heat is generally considered negligible [50]. The characteristic time constant of the perfusion heat transfer is on the order of seconds, and therefore it is more relevant for seconds-long laser exposures [51] rather than ms-short pulsed therapies. For in-vivo measurements, we can ignore potential cellular-scale variations (${\sim} $ 30 μm) in the absorption coefficient and laser-intensity variations (${\sim} $ 10 μm) in the heat source profile. since the variations would be homogenized within 1 ms of the laser pulse, as heat diffuses by about 30 μm over that duration.

5. Conclusion

We present a methodology for determination of the optical and thermal parameters in multi-layered materials, applicable to retinal photothermal therapies. We demonstrated that, by matching the optical path length changes, recorded with line-field phase-resolved optical coherence tomography, across the beam width and along the full course of heating and cooling, one can precisely identify the material properties. Based on model fittings, we estimated the temperature with a precision of about 0.2°C for a heating amplitude not exceeding 8°C during 1 ms—well within the safety range of the retina [5]. In-situ measurements of the retinal heating in every patient should allow precise titration of the energy deposition for photocoagulation and non-damaging laser therapy. We finally note that the present methodology is also generally applicable to any planar multi-layered assembly, such as skin, for example, where the method could complement the current temperature monitoring techniques, typically based on thermal cameras [52,53]. Beyond biomaterials and tissues, the method could be extended, for example, to coating and surface metrology [54] or artwork studies [55].