Interferometric thermometry of ocular tissues for retinal laser therapy

: Controlling the tissue temperature rise during retinal laser therapy is highly desirable for predictable and reproducible outcomes of the procedure, especially with non-damaging settings. In this work, we demonstrate a method for determining the optical absorption, the thermal conductivity, and the thermal expansion coefficients of RPE and choroid using phase-resolved optical coherence tomography (pOCT). These parameters are extracted from the measured changes in the optical path length ( ∆ OPL ) using an axisymmetric thermo-mechanical model. This allows the calculation of the temperature rise during hyperthermia, which was further validated by imaging the temperature-sensitive fluorescence at the same location. We demonstrate that, with a temperature uncertainty of ± 0.9° C and a peak heating of about 17° C following a laser pulse of 20 ms, this methodology is expected to be safe and sufficiently precise for calibration of the non-damaging retinal laser therapy. The method is directly translatable to in-vivo studies, where we expect a similar precision.


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][2][3][4][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][18][19][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][25][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.

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

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.

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) 1 and (r,z) 2 , at the instant t by: ( From the unwrapped phase difference, the optical path length is simply derived by: where λ OCT is the central wavelength of the OCT source.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 ∆OPL A , ∆OPL B , and ∆OPL C by: where n 0 is the index of refraction of water at room temperature T RT , u z is the vertical displacement of the plane.∆n is the average change of refractive index associated with the temperature elevation θ = T − T RT , where ∆n(θ) = n(T) − n 0 , between the planes of interest and a virtual undisturbed plane V (no temperature elevation) in the water layer located above the planes A, B, and C: See Supplement 1 for more details on derivation of the OPL changes.
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 ∆φ 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 temperaturesensitive dye can be expressed using the following model [33]: ln where T is the absolute temperature, T ref a reference temperature point of choice, and the coefficient α was found empirically (see Supplemental Document 1).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(θ(t = 20 ms)).Therefore, the temperature rise difference was calculated relative to the temperature at 20 ms, when the temperature reached its maximum value.

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]: where σ r , σ z , σ θ , are the normal stress components in the cylindrical coordinate system and τ rz is the shear stress in the r-z plane.The constitutive equations for an isotropic thermo-elastic medium are: where is the shear modulus, E is the Young's modulus, ν is the Poisson's ratio; λ = 2Gν/(1 − 2ν) is the Lamé's first parameter, ε v = ∂u r /∂r + u r /r + ∂u z /∂z is the volumetric strain.The thermal stress contribution is represented βθ, where β = 2Gα TE (1 + ν)/(1 − 2ν) is the thermo-mechanical coupling parameter, and θ is the temperature rise (θ = T − T 0 ).We note that the quasi-static assumption is not required to pursue the calculations of the tissue response.
We show in the Supplemental Document 1 that changing Eqs (9 − 10) to the wave equation does not impact the calculated ∆OPL.Therefore, we considered quasi-static conditions for simplicity.The heat flux is described by the Fourier's heat conduction law: where q is the heat flux vector q = [q r , q z ] T , κ is the coefficient of thermal conductivity and ∇ is the gradient operator.The heat flow in the z direction integrated over the time t is: and the heat diffusion equation can then be expressed as: where ρ is the material density and c p is the material specific heat capacity.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]: where P 0 is the incident laser power, a is the laser radius, ϕ laser (r, t) is the spatio-temporal profile of the laser pulse, µ a is the optical absorption coefficient, µ s ′ = µ s (1 − g) is the reduced scattering coefficient, with µ s being the scattering coefficient and g the anisotropy factor.The optical absorption coefficient is proportional to the absorber concentration c (i.e., melanin concentration) with µ a = εc, where ε is the molar attenuation coefficient.
The ϕ laser (r, t) can be decomposed into a spatial profile, i.e., the image of the fiber output, and a temporal profile: We measured the spatial profile by imaging the beam at the sample plane.The profile can be described by a convolution between a Gaussian blur (equivalent to a defocus) and a Heaviside function (disk): where H is the Heaviside function and σ blur is the Gaussian blur [36].a laser and σ blur were obtained by fitting the measured profile to the above function (see Supplemental Document 1).The temporal profile of the pulse is described by a product of a double exponential function and a Heaviside function.This expression was found to best fit the temporal profile of the pulse measured with a photodetector.
where A 1 , A 2 , t 1 , and t 2 were determined by fitting the measured pulse profile (see Supplemental Document 1).
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: q e (r, t) = P 0 πa 2 ϕ laser (r, t) where l is the thickness of the absorbing layer.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: Q e (r, t) = ∫ t 0 q e (r, τ) dτ. (23)

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 [ũ 1 r , ũ0 z , θ0 ] T , where ũ1 r and ũ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 θ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 θ(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 θ(r, z, t), the local change in the index of refraction ∆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, ∆n is calculated by averaging ∆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 ∆n and u z (r, z, t) calculated, the OPL changes are fully determined (Eqs (4 − 6)).

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]: where A is the rate constant, E* is the activation energy, R is the gas constant (8.314J.K −1 .mol−1 ) and T(t) is the absolute temperature as a function of time.The total decrease of N from its initial value N 0 over the course of hyperthermia with duration t h is described by the Arrhenius integral Ω: ln A (1.6 × 10 55 s −1 ) and E* (340 kJ/mol) have been calibrated in previous studies to set Ω so that a value of 1 corresponds to the RPE damage threshold [9].It has also been shown that under Ω<0.1 there was no detectable cellular response, which thus defines the non-damaging therapeutic window of 0.1<Ω<1 [15].The Arrhenius integral can be directly calculated using the temperature course at any chosen point in the tissue.

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, µ a = ε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 0 .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 ∆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.

Temperature calculations and precision
After determination of the unknown materials (α TE , κ, and c/c 0 ) from the ∆OPL data, the temperature evolution can be calculated at any point in space and time using our thermomechanical 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: Ω = 6 × 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.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).

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 ∆θ(t) = θ − θ ref according to Eq. (8).To compare the temperature measured by fluorescence with the calculated temperature profile θ calc , θ ref was set to match θ calc at 20 ms.∆θ obtained from fluorescence was hence converted to the temperature rise θ.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.

Discussion
From the temperature time course and the corresponding uncertainties, it is possible to compute the uncertainty in the Arrhenius integral Ω.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 Ω 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 Ω 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 Ω(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.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][52][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. 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.

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.

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

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

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

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

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 0 : relative melanin concentration, with 95% confidence bounds. a
www.engineeringtoolbox.com.Data taken at 25°C.References are indicated next to values. a