1. Introduction

Many laser-based treatments rely on the transient hyperthermia of tissues to trigger a therapeutic response, where the temperature increase results from the absorption of the laser radiation and conversion into heat [1–5]. In retinal treatments [6,7], visible or near-infrared light transmitted through the anterior segment and vitreous is absorbed primarily by melanosomes, which are highly concentrated in the retinal pigmented epithelium (RPE) [8] and in the pigmented choroid [9]. Therefore, the extent of tissue heating highly depends on the optical absorption coefficient, which is proportional to the absorber concentration. While the other thermal properties of ocular tissues, such as thermal conductivity and heat capacity, are typically regarded as constant across a population and as regionally-uniform within the eye, the melanin concentration can vary widely between individuals and locally in the eye by up to a factor of four [10]. Such pigmentation variations lead to large variations in the temperature rise during laser therapy. Moreover, ocular transparency decreases with age, potentially by tens of percent [11,12], which decreases the temperature accordingly.

The therapeutic response to hyperthermia is commonly characterized using the Arrhenius integral Ω, which assumes first-order reaction kinetics for the rate-limiting chemical reaction of thermal denaturation. It can be calculated by considering the temperature time course [13], where a value of Ω = 1 corresponds to the threshold of tissue damage [9]. Since it was also shown that Ω < 0.1 leads to no therapeutic effects [5], these two bounds define the relatively narrow therapeutic window of 0.1 <Ω < 1 of non-damaging laser retinal therapy [5,14,15]. For example, for an incident laser power of 13 mW, focused to a spot of 120 μm in diameter, which are typical settings in non-damaging retinal laser therapies, a 20% variation in the optical absorption coefficient can bring the treatment out of the therapeutic window (see Supplemental Document 1). It is therefore critically important to adjust the laser dosage (i.e., the laser power or pulse duration) according to the patient’s optical properties to ensure successful and non-damaging treatment.

Currently, in clinical applications, the laser damage threshold is typically measured first by visually monitoring the test lesions in the periphery of the retina [16], but this method is damaging and inherently subjective. Alternatively, opto-acoustic methods are being developed as a more objective approach to patient-specific laser power titration in the context of retinal treatment [17–20]. An acoustic transducer, placed on the cornea, detects the pressure waves generated upon the fast thermal expansion of the melanosomes under laser pulses. The technique yields a volume-averaged temperature that one needs to convert to a peak temperature using calibration [21]. Optical coherence tomography (OCT) has also been proposed as a non-invasive method to monitor retinal changes during treatment [22,23]. Beyond structural information, OCT recordings provide rich information that can give insights into tissue responses to heat [24–26]. For example, researchers have suggested using fringe washout due to bubble-induced mechanical displacement to monitor thermal effects in real-time [27]. Others have suggested using the phase information, or optical path length changes (ΔOPL), from OCT signals for a more quantitative interpretation of the thermal expansion [22,23,25,28].

The OPL changes, however, are not only related to the thermal variations of the refractive index (via the thermo-optic coefficient), but also to the mechanical displacements of the tissues associated with the thermal expansion. To interpret such data, it is necessary to relate the OCT measurements to a thermo-mechanical model that considers both these effects. This model should either consider well-determined material properties or be used to fit the unknown or poorly determined parameters. Along these lines, our group has previously proposed a method that combines measurements of ΔOPL using high-speed phase-resolved OCT signals with an analytical thermo-mechanical model [29]. The procedure was first demonstrated using a polymer-based multi-layer tissue phantom and, after parameters fitting, it yielded an expected temperature precision superior to that of other non-invasive techniques. Here, we extend this work to ocular tissues ex-vivo and validate it by a temperature-sensitive fluorescence measurement. We demonstrate that, by using the ΔOPL data, we can fit thermo-mechanical properties of the tissue, namely the optical absorption coefficient, the coefficient of thermal expansion, and the thermal conductivity. We show that, given the uncertainties in the parameter determination, the temperature precision is sufficient to ensure that laser treatments would remain within the therapeutic window. We finally discuss the applicability of the method to in-vivo retinas.

2. Materials and methods

2.1 Optical setup

The optical layout (Fig. 1) follows the design by Veysset et al. [29] and Pandiyan et al. [30] for high-speed line-field phase-resolved OCT (pOCT) imaging. A superluminescent diode (M-T-850-HP-I, Superlum, 9-mW power) with a center wavelength of 840 nm and a bandwidth of 110 nm was used for illumination, with a theoretical axial resolution of 2.8 µm in air. The beam was first collimated with a reflective collimator to a diameter of 1.73 mm. A cylindrical lens (CL1) focused the beam along the y-axis into a line at the sample plane (gaussian profile, measured FWHM = 220 µm along the x-axis), after optical conjugation by afocal telescopes L1 + L2, L3 + L4, and L5 + EML (eye-model lens). The beam was split between a reference path and a sample path after L1 using a 30:70 (R:T) beam splitter. In the sample path, a one-dimensional galvo scanning (GS) mirror (8310K Series Galvanometer, Cambridge Technology) steered the beam in the y-direction, and a deformable mirror (DFM) (DM69-15, ALPAO) was used for fine adjustment of the imaging focus. Both the GS and DFM were conjugated with the entrance pupil plane (EML plane). The OCT illumination power at the sample was 1.3 mW. The image plane was then anamorphically conjugated to the slit in the detection path using two cylindrical lenses (CL2 and CL3) [29]. Finally, the OCT spectrometer included a 600 line/mm grating (Wasatch Photonics) and a high-speed camera (Phantom v641, 10-µm pixel size) with a resolution of 512 × 768 pixels (spatial × spectral). The magnification between the sample plane to the camera was ×6.4, where 1 pixel corresponded to 1.6 µm at the sample. The full field of view of the B-scan is 800 µm laterally (x-axis) by 980 µm in depth (z-axis). B-scans were recorded at a frame rate of 10 kHz (98 µs exposure), which is the maximum camera rate for this resolution, for a total duration of 100 ms. Flat optical windows (SF11 glass, Schott) were added to the reference path to manually compensate the dispersion between the reference and the sample paths (not shown).

 

Fig. 1. Optical setup. (a) Top view (Y-Z plane) including the line-field spectral-domain OCT (red path), the heating laser (green path), and the fluorescence imaging (orange path). SLD: superluminescent diode, RC: reflective collimator, CL achromatic cylindrical lens doublet, L: achromatic lens doublet, BS: beam splitter 30:70 (R:T), M: mirror, GS: galvo-scanning mirror, DFM: deformable mirror, DCM: dichroic mirror, EML: eye-model lens. DPSS: diode-pumped solid-state laser, BPF: band pass filter, P: pupil plane, I: image plane. The effective focal lengths of the optics are: RC = 15 mm, CL1 = CL1R = 250 mm, L1 = 125 mm, L2 = 150 mm, L3 = 400 mm, L4 = 250 mm, L5 = 75 mm, EML = 10 mm, CL2 = 200 mm, CL3 = 75 mm, L6 = 250 mm, L7 = 100 mm, L1R = L2R = 250 mm, RC’ = 33 mm, L1’ = 250 mm, L2’ = 500 mm, L3’ = 300 mm, L1’’ = 300 mm, L2’’ = 100 mm. (b) Unfolded top view (Y-Z plane) of the OCT path showing the illumination beam (solid beam) and two imaging rays in the detection path (red lines). (c) Corresponding side view (X-Z plane). (b,c) Schematics showing the anamorphic imaging configuration based on cylindrical lenses (CL2-CL3).

Download Full Size | PDF

We overlapped the OCT optical path with a heating beam path using a dichroic mirror (DCM, NFD01-532-25 × 36, Semrock) such that the OCT line illumination coincided with the round heating beam spot (see Fig. 1(a)). The heating beam source was a 532-nm diode-pumped solid-state laser (85-GHS-305-042, Melles Griot) coupled into an optical fiber (NA = 0.22, diameter = 400 µm). To obtain a relatively uniform heating beam spot, the fiber output was imaged using 3 afocal telescopes (RC’+L’1, L’2 + L’3, L5 + EML) such that the object planes of the OCT and the heating paths matched. An iris was introduced at an intermediate image plane to control the size of the heating spot to yield a diameter of about 120 µm. The pulse duration and timing were electronically controlled using a delay generator (DG535, Stanford Research Systems, Inc.), synchronized with the camera acquisition. For all experiments, the pulse duration was set to 20 ms and the heating pulse was triggered 10 ms after the camera starts acquiring. We set the power to 10.5 mW in order to operate below an Arrhenius integral value of Ω = 0.1, providing us some room for initial estimation of the absorption coefficient, while remaining below the damage threshold Ω = 1, thereby allowing us to repeat the measurements in a reversible manner and average the datasets. The temporal and spatial distribution of the heating beam were measured and fitted by analytical functions to serve as an input to the simulations (see Supplement 1).

To provide an independent validation of the tissue temperature rise θ, we combined the OCT system with a fluorescence imaging setup. We used a temperature-sensitive dye ERthermAC (BioTracker, Millipore Sigma), with a peak of absorptivity around 532 nm, matching the wavelength of the heating beam [31]. This allowed us to use the heating beam also for recording the temperature-dependent fluorescence signal while the heating laser was on. We used the same high-speed camera to image the fluorescence by bypassing the OCT optics from the DFM to the grating (orange path in Fig. 1) with a total magnification of 9.3, where 1 pixel was 1.1 µm at the sample. We switched from fluorescence to OCT imaging using two mirrors mounted on magnetic mounts (between L4 and DFM, and between the grating and L7). Fluorescence images were recorded at an acquisition rate of 2 kHz (490 µs exposure).

2.2 Tissue preparation

All experimental procedures were approved by the Stanford Administrative Panel on Laboratory Animal Care and conducted in accordance with the institutional guidelines and conformed to the Statement for the Use of Animals in Ophthalmic and Vision research of the Association for Research in Vision and Ophthalmology (ARVO). Long Evans rats (n = 5) were used in these experiments. Animals were euthanized using an intracardiac injection of phenytoin/pentobarbital (Euthanasia Solution; VetOne, Boise, IDA, USA). The eyes were enucleated and rinsed in phosphate buffered saline (PBS), anterior segment and lens were removed, the eye cup was cut to a 3 × 3 mm square centered on optic nerve, and then further split into 4 leaves to allow better flattening. Retinal pigment epithelium with the underlying choroid and sclera were dissected from the rest of the eye.

The fresh tissue was then sandwiched between two glass coverslips of 200-μm and 1-mm thickness, facing the beam and the back side, respectively, and immersed in PBS. The assembly was then placed in the optical setup in an inverted microscope configuration. After the pOCT recording, the tissue was placed for 30 min at 37°C in a scintillation vial containing the temperature-sensitive dye ERthermAC dissolved in DMSO at a concentration of 5 μM, following the procedure from Brinkmann et al. [32]. The tissue was then washed with a PBS solution at room temperature and was placed again between the two glass substrates for fluorescence measurements. The OCT and the fluorescence measurements were taken at approximately the same location on the tissue, so we could assume that pigmentation and the material properties of the tissue were the same for both measurements. The sample assembly therefore consisted of a glass as the first layer facing the laser beams, a PBS layer, the RPE, the choroid, the sclera, a second PBS layer, and a second glass layer. For modeling purposes, the front glass and the second PBS layers were considered as a semi-infinite media, while the second glass layer was ignored (see Supplemental Document 1).

The study was performed with five animals. We show the full data analysis and data fitting for the Animal 1, where we determine the coefficient of linear thermal expansion α_TE, thermal conductivity κ, and the melanin concentration relative to their initial literature values. The other animals’ data was used to fit for their respective melanin concentration using the α_TE and κ values established with Animal 1.

2.3 OCT and fluorescence signals recording and processing

At each location in the explant, we averaged five heating pulses recorded by OCT and by fluorescence. For phase-sensitive OCT reconstruction [30], we followed the standard initial processing steps, including the background subtraction, wavenumber resampling, numerical dispersion compensation, and time referencing to cancel the local arbitrary phase [30]. We did not register the acquired B-scans because the sample was relatively static, as confirmed by the phase stability observed in our OCT analysis. The OCT signal processing steps listed below follow the procedures described in [29].

Since the heat source is a axisymmetric (approximately a disk), the thermo-mechanical problem is more conveniently described in polar coordinates (r,z). The OCT imaging probes the tissue response along the heating disk diameter, so the radial r-axis matches the x-axis shown in Fig. 1(a). After initial processing, which eliminated the local starting arbitrary phase, the OCT signal can be described as a complex number:

From the signal, we can compute the phase difference between two points at desired spatial coordinates, (r,z)₁ and (r,z)₂, at the instant t by:

From the unwrapped phase difference, the optical path length is simply derived by:

The amplitude of the OCT signal was used to generate B-scan images (see Fig. 2(a)) of the explant, from which we identified the RPE (z = 0) as the brightest layer (plane A). From this plane, we defined two planes, B and C, 10 and 20 μm deeper, respectively. For each plane, we used a phase reference point that was at the same depth (same z-coordinate) but radially away from the heating zone (r_ref = 150 μm), where the tissue was little disturbed, while the OCT signal still maintained a sufficient SNR to limit the phase noise. We can thus define $\Delta OP{L_A}$, $\Delta OP{L_B}$, and $\Delta OP{L_C}$ by:

See Supplement 1 for more details on derivation of the OPL changes.

 

Fig. 2. (a) B-scan of the tissue sample showing the RPE, choroid, and sclera. The planes A (z = 0), B (z = 10 μm), and C (z = 20 μm), at which we extracted the ΔOPL data, are indicated by horizontal arrows. (b). Time course of the experimental data and fitted ΔOPL at the center of the heating beam (r = 0). (c) Radial cross-section of the experimental and fitted ΔOPL data at two time points during the heating phase: t = 10 ms (left) and t = 20 ms (right). (d,e,f) Full spatio-temporal representation of the experimental (top) and fitted (bottom) ΔOPL data for the planes A, B, and C, respectively.

Download Full Size | PDF

For each point along the individual planes, the complex value of the OCT signal was averaged over 3 pixels in depth (± 1 pixel in z) and 5 pixels radially (± 2 pixels in r). The complex value of the reference point at r_ref was averaged radially over 7 pixels (± 3 pixels in r) to compensate for the lower OCT SNR at this location. For all time points, we computed, using Eq. (2), the phase difference data $\Delta \varphi $ between each point along the plane within the radius r < 120 μm and the fixed reference point on this plane (see Supplemental Document 1 for further details on data processing). We also note that it is possible to extract the OPL data at additional planes (not limited to 3) but at the expense of larger datasets.

The temperature dependence of the fluorescence intensity I(T) emitted by the temperature-sensitive dye can be expressed using the following model [33]:

Similar to OCT imaging, we triggered five laser pulses and recorded fluorescence images, where the 532-nm beam acted both as the heating beam and as the fluorescence excitation. This way, we measured changes in fluorescence intensity for the duration of the heating. The fluorescence intensity of the captured images was, firstly, normalized by the heating beam intensity, as fitted to the measured temporal profile (see Supplemental Document 1) and, secondly, averaged in each frame over an area with a radius of 20 µm from the center of the heating spot, where the temperature was relatively uniform, to obtain the temporal variation of the fluorescence intensity during heating. Since the laser power gradually increases at the beginning of the pulse and plateaus after a few ms, we picked the reference intensity point ${I_{ref}}$ to be the fluorescence intensity at the end of the pulse t = 20 ms, with ${I_{ref}} = I({\theta (t = \textrm{20 ms})} )$. Therefore, the temperature rise difference was calculated relative to the temperature at 20 ms, when the temperature reached its maximum value.

2.4 Governing equations of the axisymmetric thermo-mechanical problem

As shown in Eqs (4 − 6), to model the ΔOPL, the vertical displacements and temperature fields should be calculated and these quantities are correlated during thermal expansion. Thermal expansion can essentially be treated as a thermo-mechanical problem, with coupled equations of elasticity and heat diffusion. The governing equations we reproduce here were described in our previous paper and applied to a tissue phantom [29]. Here, we also assumed isotropic and uniform material properties within each layer and linear material behavior.

Assuming quasi-static conditions. the governing partial differential equations of equilibrium for an elastic medium in cylindrical coordinates are [34]:

The heat flux is described by the Fourier’s heat conduction law:

When laser light is absorbed in tissues, its energy is converted into heat, which can be described using the Beer-Lambert law as a heat source term [35]:

The ${\phi _{laser}}({r,t} )$ can be decomposed into a spatial profile, i.e., the image of the fiber output, and a temporal profile:

In reality, laser light intensity exhibits a continuous exponential decrease with the tissue depth, but for modeling purposes we discretized the volumetric heat sources into discrete surfacic heat sources in the radial plane (see Supplemental Document 1). For the 1-μm thick melanosome-rich layer within the RPE, the amplitude of the surface heat source was set to equal the total heat deposited over the 1-μm thickness. The surfacic heat source approximation is justified since the surface heat source and volumetric heat source become equivalent after the characteristic time of diffusion (0.7 μs) over a layer thickness of 1 μm. The surface heat source therefore simplifies as:

In a similar fashion, the absorbing pigmented choroid was discretized in our model into five 5-μm layers and a surface heat source was applied at the top of each layer, taking into account the light transmission through upper layers (see Supplement 1). This brings the number of heat source planes to a total of six. The external heat flow source at these planes is obtained by integrating the heat flux over time from 0 to t:

2.5 Solution to the axisymmetric thermo-mechanical problem

In our previous work we showed that it is possible to analytically solve the thermal expansion problem described above in the Hankel-Laplace (HL) domain. The solution is fully derived in the Supplemental Document 1. In brief, taking advantage of the planar layered structure, it consists in transforming the equations presented above into the HL domain and solve for the state vector ${[{{{\tilde{u}}^1}_r,{{\tilde{u}}^0}_z,{{\tilde{\theta }}^0}} ]^T}$, where ${\tilde{u}^1}_r$ and ${\tilde{u}^0}_z$ 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. The relationship between states vectors at different depths is found using a transfer matrix [37,38] or a stiffness matrix [39]. The numerical inverse Hanke-Laplace transformations then convert the solutions into the space and time domains to obtain the variables of interest: ${u_z}(r,z,t)$ and $\theta (r,z,t)$. This approach has been commonly applied to geophysics problems [40,41], was shown to be more than 10× computationally efficient [42,43] compared to finite-element methods (FEM), and allows the definition of semi-infinite spaces. We have validated this approach in our previous work [29] and have adapted it here to model the tissue response.

After finding $\theta (r,z,t)$, the local change in the index of refraction $\Delta n$ is calculated using the empirical form for the dependence of refractive index of water on temperature $n(T )$ [44]. For simplicity, we assume that the refractive index of the tissues follows the same temperature dependence. Finally, $\overline {\Delta n} $ is calculated by averaging $\Delta n$ computed at 40 planes between the planes of interest A, B, and C and a virtual plane V located 150 μm above the RPE layer in water (Eq. (7)), where the temperature elevation is negligible. With $\overline {\Delta n} $ and ${u_z}(r,z,t)$ calculated, the OPL changes are fully determined (Eqs (4 − 6)).

2.4 Arrhenius integral for determination of the therapeutic window

The Arrhenius model is commonly used to assess the thermal damage in retinal cells during transient hyperthermia. It assumes that cells are irreversibly damaged when thermal denaturation of a critical molecular component, which follows first-order kinetics, exceeds a certain fraction. Under this approximation, the concentration N of the critical component decreases at the rate [9]:

3. Results

3.1 OCT imaging and thermo-mechanical model fitting

The ΔOPL data for the A, B, C planes, obtained after the data filtering procedure described in section 2.3, are shown in Fig. 2(d)–(f). The ΔOPL data at the center of the beam (r = 0) as a function of time is shown in Fig. 2(b), and along the radial coordinate at two time points (t = 10 and 20 ms) in Fig. 2(c). Experimentally, we observe that for all three layers, ΔOPL decreases while the laser is on (from 0 to 20 ms). As the temperature increases, tissues thermally expand; the RPE layer, where most of the heat is deposited, moves up (${u_z} < 0$) and the refractive index of the water above the RPE decreases as it warms up. Both of these effects reduce the optical path length above the RPE (Eqs (4 − 6)). When the laser is turned off, heat diffuses away and the sample cools down to finally return to its initial state. The time courses of the ΔOPL traces are similar for all three planes but they show different amplitude, which is indicative of the thermal gradient within the tissue. The heat diffusion can also be seen spatially, where the laser profile of heat deposition (Figure S2) evolves into a Gaussian-like profile of ΔOPL (Fig. 2(c)).

In fitting the experimental ΔOPL data to the model, three parameters were adjusted: the coefficient of linear thermal expansion α_TE, the thermal conductivity κ, and the melanin concentration relative to their initial literature values. The coefficient of the linear thermal expansion has not been well characterized for retinal tissues. It is often assumed to be that of water [45,46], but recent investigations by Müller et al. suggest a higher value, up to 4 times that of water [23]. Likewise, there are variations in the thermal conductivity values for retinal tissues in the literature [47,48]. As stated earlier, the melanin concentration can vary widely between individuals and therefore should be fitted. Since the optical absorption is proportional to the melanin concentration, ${\mu _a} = \varepsilon c$, we started with the optical absorption coefficients for the RPE, choroid and sclera from Hammer et al. [49] and defined a default value for the melanin concentration c₀. We could then fit the melanin concentration c relative to this default value, assuming that relative change of melanin concentration between individuals affects all layers equally. All other thermo-mechanical properties of retinal tissues have either been experimentally measured with sufficient accuracy or are agreed upon in the community (see Supplemental Document 1 for all material properties).

The values used to compute the initial guess and the fitted values with their 95% confidence interval are listed in Table 1. First, we note a good agreement between the experimental and the fitted ΔOPL, as shown in Fig. 2(b) − f, both radially and temporally. As elaborated in our previous work [29] the fitting parameters are relatively decorrelated in the calculations of ΔOPL. The coefficient of thermal expansion impacts the first term in the ΔOPL equations (Eqs 4 − 6) through the vertical displacements but not the second term. The coefficient of thermal conductivity κ governs the heat diffusion, so it is mainly fitted through the temporal profile and the lateral extent of the ΔOPL data. The melanin concentration c influences both terms in ${u_z}$ and $\Delta n$ in Eqs (4 − 6) as it dictates the amount of heat deposited locally. Importantly, it affects the distribution of heat between the planes, through the depth, and therefore directly impacts the ΔOPL gradient from the planes A to C. Decorrelation between the fitting parameters is a pre-requisite for convergence of the fit to a global minimum. The fitted value for the coefficient of thermal expansion is 2.6× higher than the value for water, but not as high as predicted by Müller et al. [23]. The thermal conductivity is within the literature range and the relative melanin concentration is reasonably close to the default value.

Table 1. Material parameters of ocular tissues from literature and obtained by ΔOPL fitting. α_TE: coefficient of thermal expansion, κ: coefficient of thermal conductivity, c/c₀: relative melanin concentration, with 95% confidence bounds.^a

View Table | View all tables in this article

3.2 Temperature calculations and precision

After determination of the unknown materials (α_TE, κ, and c/c₀) from the ΔOPL data, the temperature evolution can be calculated at any point in space and time using our thermo-mechanical model. Since we are primarily interested in the RPE temperature during laser treatment, where the temperature is the highest and where the Arrhenius integral is evaluated, we computed the temperature profile for that layer (z = 0) (Fig. 3(a)). Figure 3(b) shows the time course of the temperature elevation at the RPE surface (r = 0, z = 0), where the temperature rise reaches its maximum of 16.5 °C at the end of the heating pulse (t = 20 ms). From the temperature course, we verified that the corresponding Arrhenius integral, calculated using Eq. (2)5 for this point, is below the therapeutic level: $\Omega = 6 \times {10^{ - 2}}$. This confirms that with these laser settings, the tissue was heated reversibly and multiple cycles of heating could be applied to average the ΔOPL data.

 

Fig. 3. (a) Spatio-temporal distribution of the temperature at the surface of RPE (z = 0). (b) Time course of the temperature rise at the center of the heating beam (r = 0) calculated with the best-fit parameters (black line) and results from Monte-Carlo (MC) simulations (n = 40 trials) considering the parameter uncertainties (orange lines).

Download Full Size | PDF

The precision of the temperature calculation can be evaluated by taking into account the parameter uncertainties derived from the ΔOPL fitting. We ran Monte-Carlo simulations [50] using the fitted parameters and their respective uncertainty (Table 1) over 1000 iterations to compute the peak RPE temperature at t = 20 ms at the center of the beam (r = 0, z = 0). These simulations yielded a temperature uncertainty of ±0.9°C. The computation of these 1000 iterations at a single point in time and space took 50 min. For comparison, such calculations with COMSOL would have taken >10 days. We also performed Monte-Carlo simulations of the whole time-course of temperature at that point in space over 40 iterations. The results illustrate the uncertainty on the temperature profile, as shown in Fig. 3(b).

3.3 Validation of the temperature calculations by fluorescence imaging

To validate the interferometric measurements of temperature and our model in general, we performed an independent measurement of the temperature profile using fluorescence imaging. As described in the Materials and Methods section, we recorded fluorescence images during the heating phase using the 532-nm beam for both heating and excitation of the fluorescence. The intensity of the fluorescence as a function of time $I(t)$ was translated to a temperature difference $\Delta \theta (t) = \theta - {\theta _{ref}}$ according to Eq. (8). To compare the temperature measured by fluorescence with the calculated temperature profile ${\theta _{calc}}$, ${\theta _{ref}}$ was set to match ${\theta _{calc}}$ at 20 ms. $\Delta \theta$ obtained from fluorescence was hence converted to the temperature rise $\theta $. Since the fluorescence intensity in these measurements is averaged over a radius of 20 µm, we calculated and averaged the temperature over the same radius. Moreover, we assumed that the dye penetrated the RPE cells and was distributed uniformly over the RPE thickness (4 µm). Therefore, we also averaged the calculated temperature over a 4-µm depth with an exponential weighting function to represent the attenuation of the fluorescence excitation and emission by the tissue. We note, however, that depth averaging has little influence on the obtained temperature since the temperature rapidly (<0.1ms) homogenizes over that thickness. The comparison between the temperatures obtained by fluorescence and from calculations based on ΔOPL is shown in Fig. 4. The good agreement between the two curves validates our model and the fitting results.

 

Fig. 4. Comparison of the RPE temperature measured by fluorescence and obtained from pOCT. The orange shaded area depicts the standard error from the fluorescence measurements.

Download Full Size | PDF

4. Discussion

From the temperature time course and the corresponding uncertainties, it is possible to compute the uncertainty in the Arrhenius integral $\Omega $. While the impact of the temperature uncertainty on the Arrhenius integral at the low-power calibration point is not of particular interest since it is deep in the non-damaging range, it is important to predict what the uncertainty of $\Omega $ would be at therapeutic laser settings. Using the fitting results from the calibration point, we can estimate the laser power that would bring the treatment within the therapeutic window. The relationship between $\Omega $ and the laser power P follows an exponential form (Eq. (2–5), assuming a linear regime, where the temperature rise is proportional to the laser power (see $\Omega (P)$ in Fig. 5). In the present case, increasing the laser power to 13.6 mW brings the Arrhenius integral to the middle of the therapeutic window. From this point (i.e., a hypothetical treatment point), we performed MC simulations (n = 40) similarly to those in section 3.2 to simulate the temperature profile of the RPE at the center of the beam, varying the parameters from their best-fit value within their respective uncertainty. The temperature profiles were then converted into the Arrhenius integral values. Statistical distribution of the Arrhenius integral is shown as a boxplot in Fig. 5. It demonstrates that the precision provided by our method satisfies the therapeutic requirements since the distribution remains within the therapeutic window bounds.

 

Fig. 5. Arrhenius integral as a function of laser power. The calibration point, used to determine the coefficient of thermal expansion, thermal conductivity, and melanin concentration, is outside of the therapeutic window. The boxplot illustrates the uncertainty in the Arrhenius integral at the laser power corresponding to the middle of the therapeutic window.

Download Full Size | PDF

Application of our method in-vivo requires a few modifications. Contrary to the present sample, one would need to add the neural retina to the multi-layer model. Adding layers to the model is trivial in terms of computation and would only lead to longer computations. Direct measurements of the layer thicknesses by OCT also serve as inputs to the model, and have been shown to provide sufficient SNR to ensure image registration and yield high phase stability and sensitivity in patients [51–53]. As a first approximation, as shown here, we can assume a uniform coefficient of thermal expansion and a uniform thermal conductivity for all retinal layers. If we reasonably assumed that the values that we determined in this work are also valid in vivo, it would eliminate two fitting parameters that were not previously determined in the literature, α_TE and κ. The relative melanin concentration will obviously remain a fitting parameter. Following this logic, we fitted the melanin concentration in five specimens from different animals, assuming all other material properties known. The results are shown in Table 2 and the fitted ΔOPL curves can be found in the Supplemental Document 1. The obtained melanin concentrations are relatively close, given that the animals are from the same line and about the same age. We can, however, expect more variations in retinal pigmentation in humans under clinical conditions, as explained earlier [10]. It is important to point out that precision listed in Table 2 is higher than with the Animal 1 (Table 1) because the fitting procedure for Animals 2-5 assumed fixed values for all other material properties, unlike the case with Animal 1.

Table 2. Fitted relative melanin concentrations (c/c₀) with 95% confidence bounds for different animals using the coefficients of thermal expansion and thermal conductivity determined from Animal 1.

View Table | View all tables in this article

In-vivo measurements will involve additional unknows. First, the light transmission from the cornea to RPE is generally unknown and can vary between patients. We expect that this can be fitted by matching the experimental ΔOPL to the model independently from the melanin concentration c because c impacts not only the heating of RPE but also the gradient of heat distribution from the RPE to choroid. Monitoring the ΔOPL of additional layers, available in vivo, will decorrelate the fitting of the ocular transmittance from the melanin concentration.

Second, the angular magnification of the eye affects the spatial extent of the heating beam and consequently the incident laser fluence at the RPE. In our previous work with a tissue phantom [29], we have shown that this quantity can also be considered a variable parameter, along with defocus, and can be fitted through analysis of the ΔOPL. Adding these parameters to the fitting procedure would certainly lead to a larger uncertainty in the temperature calculations, but this could be compensated by further averaging and by including additional ΔOPL data extracted from additional planes. With more data to fit, the computational time to convergence could also increase, although this issue could be resolved by implementing parallel computing strategies. One could also imagine the calibrations measurements performed by pOCT in patients during their diagnostic imaging session, and then being used for treatment planning a few days later. It should finally be noted that the method is not limited to line-field OCT but is also applicable to point-scan or full-field OCT methods, which also allow phase-sensitive measurements.

5. Conclusion

We have demonstrated that from the changes in the optical path length during the thermal expansion observed with pOCT, one can determine the coefficients of thermal expansion and thermal conductivity and the melanin concentration in ocular tissues. This method was further validated by an independent measurement using a temperature-sensitive fluorescent dye. With the temperature uncertainly under 1°C and the peak temperature rise below 17 °C, this methodology is sufficiently precise for calibration of the non-damaging retinal laser therapy. The method is directly translatable to in-vivo applications, where we expect to achieve similar precision. It should increase reproducibility and safety of the retinal laser treatments thanks to patient-specific titration.