DEPTH-RESOLVED ANALYTICAL MODEL AND CORRECTION ALGORITHM FOR PHOTOTHERMAL OPTICAL COHERENCE TOMOGRAPHY By Maryse Lapierre-Landry Thesis Submitted to the Faculty of the School of Engineering of Vanderbilt University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE in Bio

Photothermal OCT (PT-OCT) is an emerging molecular imaging technique that occupies a spatial imaging regime between microscopy and whole body imaging. PT-OCT would benefit from a theoretical model to optimize imaging parameters and test image processing algorithms. We propose the first analytical PT-OCT model to replicate an experimental A-scan in homogeneous and layered samples. We also propose the PT-CLEAN algorithm to reduce phase-accumulation and shadowing, two artifacts found in PT-OCT images, and demonstrate it on phantoms and in vivo mouse tumors.


Introduction to pre-clinical molecular imaging
Molecular imaging is an essential tool to visualize the structure but also the function of biological components [1,2].It is used both in pre-clinical research on animal models and in clinical practice [3].Endogenous and exogenous molecular contrast agents can be used and imaging can be performed at all resolution scales to identify biomarkers, visualize biological interactions, or highlight a cellular pathway.In particular, researchers in the area of drug delivery need in vivo molecular imaging modalities with a wide array of contrast agents, different penetration depths and different resolution scales to study how the drug molecules or carriers distribute in the body, the organ or the cell [4,5].There is still a need for a low-cost, noninvasive, three-dimensional molecular imaging technique appropriate for small animal imaging at the organ level and allowing multimodal imaging.Current pre-clinical molecular imaging techniques are described here to highlight this need.
This technique is widely used in pre-clinical research and in the clinic.However, this method relies on expensive instrumentation and requires high concentration of exogenous contrast agents for molecular imaging [1].In comparison, CT scanners are also widely available in laboratories and in the clinic and are relatively cost-efficient.CT uses X-rays to probe the sample and images are reconstructed in 3D in post-processing.Similarly to MRI, CT also suffers from poor sensitivity to molecular contrast agents and high concentrations are thus used.Finally, both PET and SPECT use radionuclides to create contrast in the sample and a wide array of contrast agents can be used in combination with those two methods.However, cost and availability are an issue for both technologies, and even though they are used in the clinic, factors such as the presence of radiation or the need for a particle accelerator impede on PET and SPECT usage for pre-clinical research.Additionally CT, PET and SPECT all require the use of ionizing radiation.
In between whole body techniques and microscopy, ultrasound based techniques occupy an intermediate position with superior resolution (<100 μm) to whole-body techniques but inferior penetration depth (~10 cm).Basic ultrasound imaging [17][18][19] uses a transducer to send and receive sound waves (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) in the body.The sound waves are reflected by different structures in the sample and by measuring how long it takes the signal to return to the transducer the distance from the transducer to the sample structures can be calculated.Ultrasound transducers are commercially available and low-cost, which facilitate their uses in pre-clinical research and in the clinic.However, few molecular contrast agents (microbubbles) are available, which severely limits the uses of this technology for molecular imaging [1].Another ultrasound based method, photoacoustic imaging [20,21], uses a pulsed laser to heat absorbers present in the sample, and through the photoacoustic effect transforms the local thermal expansion into acoustic waves.The acoustic waves are then detected by an ultrasound transducer.Photoacoustic images are of similar resolution and can be acquired at similar depths as ultrasound imaging.For the photoacoustic effect to take place, a contrast agent that absorbs light and transforms it into heat is necessary.Endogenous agents, such as blood [22] and melanin [23], can be visualized with photoacoustic imaging, but a large array of functionalized exogenous contrast agents can also be used for molecular imaging.At this time photoacoustic imaging is only used in preclinical research.One of the limitations of both ultrasound-based methods is the requirement for the ultrasound transducer to be in contact with the sample, which can complicate the imaging procedure in itself, but also makes multimodal imaging difficult.This is an obstacle for eye imaging, where current researchers are unlikely to adopt an ultrasound-based method.
Optical methods are also important in pre-clinical molecular imaging as they are noninvasive, without ionizing radiations and inexpensive in comparison with whole body techniques.A large array of fluorophores and similar contrast agents has also been developed for those methods.Fluorescence microscopy is the most commonly used technique to achieve high resolution images, with a large selection of exogenous and endogenous contrast agents [24,25].
However, it is limited by its imaging depth (< 200 μm) and field-of-view (requires mosaic past ~10mm 2 ), and it is most commonly used in vitro and ex vivo.On a different scale, bioluminescence imaging [26,27] can be used on large field-of-views (>30 cm 2 ) with low resolution and traditionally no depth information (new systems now allow 3D imaging [28]).However, most in vivo applications of bioluminescence require animals with genetically encoded luciferase, which limits the use of this imaging technology in the clinic.

Contrast-enhanced optical coherence tomography
There is still a need for a low-cost, non-invasive, three-dimensional molecular imaging technique appropriate for small animals at the organ and intra-organ level and allowing multimodal imaging.The affordable cost is important for the technology to be accessible at the pre-clinical level and in the clinic.A non-invasive imaging technique allows for longitudinal studies on the same animal cohort and cause minimal disturbance to the model being imaged.
Both of those requirements will contribute to bridging the gap currently seen between results obtained in vitro and ex vivo and results obtained in human patients.Three dimensional imaging with a large field-of-view and high resolution provides an opportunity to study heterogeneous distribution at the organ level.Researchers in the area of drug delivery, for example, aim for the homogeneous distribution of their new drug or drug carrier inside a target organ or tissue.Three dimensional imaging is also important for retinal imaging, since the retina is composed of many distinct layers across its thickness.Finally, multimodal imaging creates a more complete picture of the sample by allowing simultaneous imaging of tissue structure, vasculature and the molecular target.This is particularly important in drug delivery studies where the target of interest is delivered to the tissue through the vasculature, and all three components (tissue, blood vessel, drug) can be imaged at once.A technology based on optical coherence tomography (OCT) could fulfill those requirements if an appropriate contrast mechanism is identified OCT is an optical method that uses near-infrared low coherence light to image tissue structure in three dimensions [29,30].Using a similar concept to ultrasound, a pulse of light is sent into the sample and backscattered at a certain location in the tissue before going back to the detector.Each pulse of light is split in between two paths of equal length: the reference arm containing a mirror and the sample arm containing the sample that is being imaged.
Interferometry is used to calculate the difference in path between photons that were backscattered in the sample and photons that went down a reference arm with a mirror.OCT offers high contrast between tissue structures of different index of refraction.It has become the standard of care in retinal imaging in the clinic [31] and is frequently used in pre-clinical research to image the eye [32], the skin [33], the brain [34], blood vessel walls [35], and tumor xenografts [36].In addition to tissue structure, OCT can provide information on blood vessel location [37], blood velocity [38] and collagen structure [39].However, OCT in its traditional form is not a good candidate for molecular imaging as the source of contrast in OCT is scattering and change in index of refraction.Contrast agents can be used to change tissue attenuation, but their detection is non-specific and is not appropriate for high scattering environments such as tumors and lesions on the retina.Multiple techniques have been proposed to enable molecular imaging in OCT.
Spectroscopic OCT (SOCT) [40] detects the presence of a contrast agent by analyzing the absorption spectra of the tissue as it is illuminated by the broadband OCT source.Contrast agents such as melanin [41], Hb/HbO 2 [42], some commercially available near infra-red dyes [41], gold nanorods [43] and gold nanocages [44] have all been used for SOCT, which thus has the potential for molecular imaging.SOCT also has the advantage of being low-cost and to require no other equipment than a commercial OCT system.However, SOCT provides the best results when the attenuation profile of the sample is known, which can be impractical for in vivo experiments, especially longitudinal studies in animals.Additionally, the absorption peak of the contrast agent needs to be in the bandwidth of the OCT source, which will attenuate the signal used for imaging.This causes a reduction in imaging depth compared with traditional OCT.
Magneto-motive OCT (MM-OCT) [45] is an active detection method that causes a perturbation in the sample in order to detect the presence of contrast agent.Magnetic microparticles are injected in the sample and a variable magnetic field is applied to the sample while it is being imaged by an OCT system.Under the effect of the magnetic field, the particles oscillate slightly in the sample, which is then detected as a change in phase in the OCT signal.[45], live tadpoles [46] and rat mammary tumors in vivo [47].Although MM-OCT can successfully identify the contrast agent location, it is limited in its ability to quantify the contrast agent's concentration.Since the oscillation of the magnetic particles depends on the viscoelastic properties of the tissue, a higher signal could be produced both by an increase in concentration of the contrast agent or a change in the tissue's elasticity.MM-OCT is also limited in imaging speed since the strongest oscillations can usually be observed around 50Hz, and multiple oscillations needs to be acquired for each scan.

MM-OCT has been used to image cells in vitro
Pump-probe OCT (PP-OCT) [48][49][50] uses a pulsed laser to excite or bleach a fluorophore in the sample.Fluorophore with transient absorption are chosen, i.e. they are excited by the pump beam to an excited state where they absorb the OCT light.Fluorophores that have transient absorption in the near infra-red are chosen so to overlap with the OCT detector.By collecting multiple A-scans at the same position with the pump alternatively on or off, the difference in OCT signal is analyzed and the location of the contrast agent is determined.PP-OCT enables in vivo molecular imaging but requires an expensive pulsed laser, is limited to a few eligible dyes and relies on imaging short-lived excited states which are not always easy to capture.
Photothermal OCT (PT-OCT) is an additional method that has been developed to achieve molecular contrast with OCT.We aim to improve PT-OCT so it addresses the shortcomings of other contrast-enhanced OCT methods and constitutes a practical OCT-based molecular imaging technique.

Photothermal optical coherence tomography as a molecular imaging technique
PT-OCT has steadily been receiving interest since its introduction less than a decade ago [51,52].Its instrumentation is simple, requiring only an additional amplitude-modulated laser to be co-localized to the OCT imaging beam.To detect the presence of contrast agent, the amplitude-modulated laser is tuned to the absorption peak of the particle or molecule and directed toward the sample.The contrast agent will then absorb the light repeatedly and heat the sample locally by a few degrees.The resulting sample expansion and change in index of refraction will result in a change in optical path length detected by the OCT system.For most samples, the change in temperature will be less than 5°C, which lead to a change in optical path length on the order of nanometers.This change in optical path length will cause a change in phase in the complex OCT signal.As multiple heating cycles are performed by the amplitudemodulated photothermal laser, repeated OCT scans are needed at the same sample location over time to observe the cyclical change in phase in the OCT signal and determine the position of the contrast agent in post-processing.
PT-OCT has been demonstrated in vitro [53,54], ex vivo [55,56] and in vivo [57,58] on a multitude of contrast agents, such as gold nanospheres [52], gold nanoshells [59], gold nanorods [57,58], carbon nanotubes [54], indocyanine green [60,61] and blood [62,63].As more nanoparticle designs and functionalization are investigated, it is important to have a technology that can image their distribution with high resolution and in 3D in high scattering environments.PT-OCT is a low-cost option, as commercial OCT systems can be used and the photothermal laser does not need to be pulsed (the amplitude modulation is done with a mechanical chopper or an acousto-optics modulator).It also allows for multimodal imaging with most other functional OCT modalities, such as Doppler or speckle variance OCT for blood vessels visualization.

Instrumentation
A commercial spectral-domain OCT system can be altered for photothermal imaging (see Chapter 2, Figure 2 for diagram).A superluminescent diode imaging laser with λ 0 =860 nm central wavelength and 40 nm bandwidth is used as the OCT imaging beam.The light from the SLD is split between the sample and reference arms using a 50:50 fiber coupler.Photothermal heating is provided by a wavelength tunable Titanium:Sapphire laser.The wavelength of the PT laser is tuned to match the absorption peak of the contrast agent used.The PT laser intensity is modulated by an acousto-optic modulator following a square wave with frequency f 0 =500 Hz.
The light of the PT laser is combined to the superluminescent diode and sent to the reference and sample arms using the 50:50 fiber coupler.A telecentric lens is used to focus the light onto the sample (1/e 2 beam diameter ω 0 =28 μm in air).The returning light is directed to the spectrometer and CCD integrated to the Bioptigen system using a circulator (AC photonics).For a typical scan pattern, repeated A-scans are collected at the same sample position over time (M-mode scan, Δt=100 μs integration time per A-scan) before moving to another position.Scans are collected one B-scan at a time in the sample (100 point-positions per 1mm of B-scan) and multiple Bscans are collected to form volumes (C-scans).All the raw data is saved to memory before starting data processing.

Signal analysis
In spectral-domain OCT, the signal from the reference and sample arm interferes before reaching a spectrometer and a CCD that acts as a detector.The broadband light is then spatially separated based on its wavelength.For data processing, the first step is to re-sample the data so it is linear as a function of the wavenumber k instead of being linear in wavelength λ.This is accomplished with linear interpolation and a calibrated lamp spectrum.Digital dispersion compensation is then applied to the data (to correct for dispersion mismatch in between the sample and the reference arm) [64] then a fixed pattern removal is applied to remove DC signal [65].A Chirp Z transform is then applied to the OCT spectrum to transition from wavenumber to depth as the independent variable.From the result of the Chirp Z transform is also obtained the magnitude and phase of the signal.From the intensity signal, an OCT image can be reconstructed to show the sample structure.However, the phase information is needed to obtain PT-OCT data.
If contrast agent is present at a specific depth inside the sample the absorption coefficient will be much higher than for the surrounding sample and the phase signal at the same depth will vary over time at a frequency equal to the amplitude modulation frequency of the photothermal laser.The derivative of the phase signal with time δΦ is first taken to center the phase to zero and avoid phase wrapping [57].
Re ( , ) ( , 1) A Fourier transform of the phase derivative is taken to go from the time domain to the frequency domain.If contrast agent was present at this depth inside the sample, a peak is expected on the spectrum at the frequency matching the photothermal laser modulation.The height of this peak is taken as the photothermal signal.The average amplitude of the nearby (offpeak) frequencies is taken as the photothermal noise level.This process takes place at all depth and for every lateral position along a B-scan.To display the results, the photothermal noise is often subtracted from the photothermal signal.

Introduction
Molecular imaging is widely used in pre-clinical studies, for example in drug discovery [5], and to study biological processes [1].Multiple imaging modalities have been developed for in vivo molecular imaging, each on different field-of-view and resolution scales, and each with their own advantages and disadvantages.Whole body imaging techniques include magnetic resonance imaging (MRI) [8], computed tomography (CT) [10] and positron emission tomography (PET).These whole body imaging techniques have no limit on imaging depth, but are usually high cost methods, with poor sensitivity to the molecular contrast agent (CT, MRI) and/or involve the use of ionizing radiation (CT, PET).Ultrasound and photoacoustic [21] imaging have a penetration depth of up to 10 cm and have the advantage of being low cost, with equipment already available in the clinic.Ultrasound can image contrast agents such as microbubbles [66] and some nanoparticles [67] while photoacoustic imaging can image absorbers such as blood, melanin and a wide array of metallic nanoparticles.However, both ultrasound and photoacoustic imaging do not allow simultaneous multi-modality imaging since they use an ultrasound transducer, which complicates the coupling of multiple imaging techniques in one instrument and can also complicate the coupling of the instrument to the subject.Finally, optical techniques are limited to low imaging depths but allow for both large field-of-views with low resolution (e.g.whole body bioluminescence imaging [27]) and small field-of-views with very high resolution (e.g.fluorescence microscopy [24]).
There is a need for an optical molecular imaging technique with a large field-of-view and high resolution to fill the gap between microscopy and bioluminescence imaging, and complement non-optical techniques.Photothermal optical coherence tomography (PT-OCT) is an emerging technique [51,52] that would allow the imaging of contrast agents in vivo with a field-of-view (>25mm 2 ), imaging depth (>1mm) and resolution (<15μm transverse) comparable to traditional OCT.In PT-OCT, a laser is used to heat the contrast agent embedded in the sample, which in turn heats the sample locally, and thus causes an elastic expansion and a change in the index of refraction of the medium.This causes a local variation in optical path length that is then detected by a phase-sensitive OCT system.PT-OCT thus offers similar sources of contrast as photoacoustic imaging, while allowing multiplexing with other functional OCT modalities (such as speckle variance OCT [58]).Additionally this method does not require the instrument to be in contact with the sample.
PT-OCT has been demonstrated in vitro [52], ex vivo [55] and in vivo [57] on a multitude of contrast agents, such as gold nanospheres [52], gold nanoshells [55,68], gold nanorods [53,[56][57][58], carbon nanotubes [54], indocyanine green [60,61] and blood [62,63].However, PT-OCT images are still affected by two artifacts: phase-accumulation (signal becoming brighter with depth) and shadowing below high intensity regions.A method has been proposed to correct those artifacts [69], but by taking a derivative of the PT-OCT signal this method also drastically decreases the signal-to-noise ratio of the image.
A theoretical model of PT-OCT would allow for a better understanding of the photothermal signal, give indications on how to optimize the imaging procedure, and would be a tool toward solving the problem of phase accumulation and shadowing.Separate components of the photothermal process have been modeled in the past, including the sample heating through the bio-heat conduction equation [57], the change in optical path length due to the change in temperature [69], the reflectivity of a sample with specific optical characteristics [70] and the phase sensitivity as a function of the signal-to-noise ratio in OCT [71].
We propose a model that assembles those individual components and includes the evolution of the photothermal signal with depth.This is the first model to replicate reliably a PT-OCT A-scan for a homogenous sample and for a layered heterogeneous sample.This model can be customized to the specific instrument (spectral-domain or swept-source OCT), imaging parameters, and sample characteristics.Furthermore, we have used the model to implement a variation of the CLEAN algorithm [72] which successfully eliminated the effect of phase accumulation and shadowing in phantom images and in images of heterogeneous tumors acquired in vivo.This paper will present the equations necessary to implement this analytical model of PT-OCT and our version of the CLEAN algorithm.Experimental validations in phantoms will also be presented in each section.

Theory
Model of photothermal signal generation PT-OCT relies on photothermal heating of contrast agents by a laser illumination source.
For the contrast agent to be heated, the photothermal laser light must first reach the contrast agent.The photothermal laser power P(z) was generated using a Monte Carlo simulation [73] that took into the absorption and scattering coefficient of the sample, respectively μ a and μ s (cm -1 ) and z (μm), the distance travelled within the sample.Once the power that reaches the contrast agent is known (and the laser spot size is known), we can calculate how much heat will be emitted.The bio-heat conduction equation is used to model how the temperature will vary in the sample surrounding the contrast agent when it is illuminated.The bio-heat conduction equation with a heat source is the following [57]: where T (K) is the temperature, t (s) is the time, φ (W/m 2 ) is the photothermal laser fluence rate, ρ (kg/m 3 ) is the density of the sample, c (J/kg*K) is the specific heat of the sample and α (m 2 /s) is the thermal diffusivity of the sample (α = k/ρc where k (W/m*K) is the thermal conductivity of the medium).In the case where the laser spot size is small compared to the absorption depth, radial heat transfer dominates the heat conduction equation, which can then be solved in cylindrical coordinates.The temperature ΔT at the center of the laser beam (r = 0) will rise and fall over one modulation period (2t L ) of the photothermal laser in the following way [57]: where t L (s) is the dwell time of the photothermal laser on the sample before an acousto-optic modulator reduces the laser power to zero (square waves modulation, 50% duty cycle) and is the 1/e 2 beam radius (ω 0 (μm) is the waist size inside the tissue estimated with a Monte Carlo simulation [74] and z 0 (μm) is the Rayleigh range of the objective).The change in temperature will cause the sample to undergo two changes: elastic expansion and change in index of refraction.Combined, these phenomena will cause a change in the optical path length (ΔOPL) detected by the OCT system.The change in OPL over a set amount of time (Δt) at each depth z can be calculated as follows [69]: where T 0 (K) is the initial sample temperature, n is the index of refraction, dn/dT (K -1 ) is the thermo-optic coefficient that is assumed to be constant for this model and β (K -1 ) is the thermal expansion coefficient of the sample.In our sample, dn/dT>0 and β<0 thus both effects act against each other.Over multiple photothermal excitation cycles, the variation in OPL will cause a change in phase ΔΦ (mrad) in the OCT signal, which will be detected by our instrument [71]: where λ 0 (nm) is the center wavelength of the OCT laser and the average index of refraction for the sample n = n(T 0 ) is used.

Model of photothermal signal detection
Once the photothermal signal is generated inside a sample, the signal is detected using an OCT system.The signal I(z), collected from one A-scan with a sub-resolution change in reflector position ΔΦ(z), is calculated as [71]: where δ (A/W) is the detector responsivity, e (C) is the electronic charge, S (W) is the power of the OCT laser source at the sample, Δt (μs) is the time to acquire one A-scan, R s (z) is the reflectivity at the sample as a function of depth inside the sample, and R R is the reflectivity of the reference arm.The sinc function accounts for the spectrometer fall-off with depth, and the scaling constant a is determined experimentally.We use the following equation to calculate the reflectivity of the sample arm R s (z) [70]: where μ b (cm -1 ) is the backscattering coefficient of the sample, L c (μm) is the coherence length of the OCT laser and f (μm) is the focal position of the OCT beam.Reflectivity models specific for phase-sensitive detection [70,71] were chosen over other models [75], because PT-OCT is a

phase-sensitive technique
The experimental signal also contains additive, uncorrelated Gaussian white noise, divided here into shot noise A shot (z) [71] and environmental noise A env (z): () ) where φ rand and φ' rand are uniformly distributed random phases between -π and π and A 0 is an experimentally determined magnitude that includes other sources of noise (electronic, vibrational, thermal, etc.).The experimental phase change ΔΦ exp (z) is then calculated as the phase of the signal and noise terms added together.
The experimental phase change ΔΦ exp (z,t) is acquired over time (in practice 1000 repeated A-scans at the same location inside a sample), then a Fourier transform is used to go from the time domain to the frequency domain.
The photothermal signal is then defined as the magnitude of the signal at the laser modulation frequency f 0 =500Hz, The photothermal signal is also weighted by the OCT laser wavelength λ, the average index of refraction n and the laser modulation frequency f 0 in order to have units of nm (same as ΔOPL).(See [57] for full derivation).
The photothermal noise is calculated following the same equation but where p(z,f 0 ) is replaced by the mean magnitude of nearby frequencies.

Model for heterogeneous samples
The model described in the previous sections is only valid for homogenous samples.
However, quasi-heterogeneous samples can be created by superposing layers of different homogenous samples.First the power P(z) is simulated for a multi-layer sample using the Monte Carlo method [73].Then the temperature change ΔT and reflectivity R s (z) are calculated for each layer separately and combined as one continuous function vs depth.The beam waist ω(z) is again estimated using a Monte Carlo simulation of a focused Gaussian beam in a multi-layered sample [74].The following calculation steps to obtain the PT-OCT signal are then the same as in the one-layer case.

PT-CLEAN algorithm to eliminate phase accumulation and shadowing
The CLEAN algorithm was invented in the 1970s with applications in radio astronomy [76] before being adapted for microwave imaging [77].It was then adapted for OCT to reduce speckle noise, cancel coherent artifacts and improve resolution [72,78].The CLEAN algorithm relies on a basic procedure: deconstruct an image point-by-point using the point-spread function (PSF) proper to the instrument, then re-construct the image point-by-point using a modified PSF to improve image quality.In the case of PT-OCT, we want to correct for phase accumulation and shadowing, thus a PSF without those two artifacts will be used for reconstruction.
First, our analytical model of PT-OCT is used to create the PSF of our system.Using the theoretical model for this step is more practical than to image a thin, homogeneous absorber (<10 μm) that could be approximated as a point.For the PSF, a 7 μm thick absorbing and scattering layer of tissue (thickness on the order of the OCT axial resolution) is simulated and serves as a 1D PSF of PT-OCT, denoted h(z).This PSF can be generated with approximated values for the scattering and absorption coefficients since it will be normalized in the next step of the algorithm.Second, the original experimental PT-OCT B-scan is considered, and the first pixel in depth to have a photothermal signal above a certain threshold (ΔOPL = 2 nm) is selected (coordinates [x (1) , z (1) ]).Third, the following deconvolution kernel is subtracted from the original image [72]: where ε<1 is the loop gain, or fraction of the image that is removed at each iteration of the algorithm, h[z-z (1) ] is the PT-OCT PSF centered at [x (1) , z (1) ], D image [x (1) , z (1) ] is the intensity of the selected pixel in the original PT-OCT B-scan and h max is the maximum value of the PSF.A point of intensity δ[x (1) , z (1) ] = ε•D image [x (1) , z (1) ] at position [x (1) , z (1) ] is then saved for the final image reconstruction.
The next iteration then starts by finding the first pixel in depth to have a photothermal signal above threshold in the remaining image (original image minus deconvolution kernel): then we proceed as described for the previous iteration.The iterative process is stopped when the brightest pixel is below the noise floor of the original image.The final CLEANed image is obtained by convoluting each points δ[x (i) , z (i) ] with a Gaussian function of 7 μm in width at half maximum.What is left of D image after the iterative process can be added to the CLEANed image to provide a realistic-looking noise floor.The PT-CLEAN algorithm takes into account the decreasing SNR of the image with depth since the PSF also has a decreasing SNR with depth.
Additionally, a speckle-correction algorithm can be used if necessary in combination with the PT-CLEAN algorithm to improve data visualization.
A summary of the PT-CLEAN algorithm can be seen in Fig. 1.

Instrumentation
A commercial spectral-domain OCT system (Bioptigen, Inc.) was altered for photothermal imaging (Fig. 2).A superluminescent diode imaging laser (SLD, Fig. 2) with λ 0 =860 nm central wavelength and 40 nm bandwidth (Superlum) was used as the OCT imaging beam.The light from the SLD was split between the sample and reference arms using a 50:50 fiber coupler (50:50, Fig. 2).Photothermal heating was provided by a wavelength tunable Titanium:Sapphire laser (Coherent).The wavelength of the PT laser was tuned to match the absorption peak of the contrast agent used (λ PT =750 nm for indocyanine green and brown human hair, λ PT =740 nm for gold nanorods).The PT laser intensity was modulated by an acousto-optic modulator (Brimrose) (AOM, Fig. 2) following a square wave with frequency f 0 =500 Hz.The light of the PT laser was combined to the superluminescent diode and sent to the reference and sample arms using the 50:50 fiber coupler.A telecentric lens (Bioptigen) was used to focus the light onto the sample (1/e 2 beam diameter ω 0 =28 μm in air, 36 μm simulated in tissue).The returning light was directed to the spectrometer and the CCD integrated to the Bioptigen system using a circulator (AC photonics).

Phantom fabrication
Clear silicone (Quantum Silicones) was used to create solid homogeneous phantoms.The silicone consists of two parts (base and curing agent) that need to be mixed in a 10:1 ratio.For all phantoms, rutile titanium dioxide (TiO 2 , <5μm, Sigma Aldrich) was first added to the silicone base for a final concentration of 2.05mg/g (μ s ≈50cm -1 ).The base and TiO 2 were then mixed for two minutes and degassed for two minutes using a planetary centrifugal mixer (Thinky USA).In phantoms where indocyanine green (ICG) was added, it was first dissolved into 800 μL of 70% ethanol then added to the curing agent of the silicone.It was then mixed and degassed for two minutes each in the planetary centrifugal mixed (when ICG was not used, this step was skipped).
Finally, the base and curing agent were mixed together and degassed.This final mixture was then poured into a petri-dish and placed under vacuum at 29 inches of Hg for five minutes with brief returns to standard pressure every minute.
To make layered phantoms, the final silicone mixture for each of the layers was poured between two glass slides with 115 μm spacers attached on both sides.The layers were then cured at 70°C for 12h and superimposed after curing.For single-layer phantom the silicone mixture was poured into a petri-dish and then cured following the same procedure.
A phantom was made using a human brown hair as a thin absorbing layer.A silicone mixture was prepared as described above (no ICG).Then, a single hair was attached to both sides of the petri-dish so it would form the hypotenuse of a right triangle of known height (petri dish depth) and length (petri-dish diameter).The silicone mixture was then poured into the petri-dish, covering part of the hair strand.Using trigonometry, it was thus possible to know the depth of the hair as it progressively penetrated the silicone.The phantom was then cured.

Imaging parameters and signal processing for phantom experiments
For every phantom experiment, the SLD power at the sample was 2 mW, with a 36 μm 1/e 2 beam radius at the focus (in silicone).The photothermal laser power at the sample was varied in between 2 mW and 16 mW depending on the experiment.For each 1D depth-resolved PT-OCT A-scan, 1000 OCT A-scans were recorded over time with every OCT A-scan lasting Δt=100 μs.Each OCT-A scan was resampled from a linear wavelength space to a linear wavenumber space, then corrected for dispersion and background subtracted [57].A Chirp-Z transform was used to convert the wavenumber data to depth-resolved data.The phase data at each point in depth was then used for the rest of the analysis.The temporal derivative of the phase at each point was calculated and a Fourier-transform was taken in the temporal dimension.
The data was then treated exactly like the data generated by the analytical model (see Eq. 13).
For each experimental validations of the model, one B-scan (4 mm in length) containing 400 PT-OCT A-scans was averaged to produce one noise-reduced PT-OCT A-scan and facilitate data visualization .The same method was used with the data generated by the model.

In vivo imaging
All animal work was approved by the Vanderbilt IACUC committee, and all surgical procedures were performed using aseptic techniques.One nude female mouse (Foxn1 nu / Foxn1 nu , The Jackson Laboratory) underwent dorsal skinfold surgery to enable imaging of MDA-MB-231 (human, breast) tumors.Following baseline (pre-injection) PT-OCT imaging (no signal detected), gold nanorods (AuNRs) coated with methoxy-terminated poly(ethylene glycol) (PEG) were injected (200 μL injection, 9 nM) with a retro-orbital injection.The mouse was imaged again 17h post-injection with PT-OCT while under anesthesia (isofluorane mixed with air, 1.5-2%).The full protocol, including AuNRs synthesis and surgical procedure is described in detail in [58].For this experiment, the photothermal laser was set to λ PT =740 nm, with 70 mW of power at the sample.The SLD delivered 6.6 mW of power at the sample.B-scans were acquired and processed as explained in the previous section.

Input parameters for the analytical model
Two different samples were simulated using the PT-OCT theoretical model: a one layer, 5 mm thick silicone phantom containing TiO 2 (2 mg/g) and ICG (0.09 mg/g), and a silicone phantom containing TiO 2 (2 mg/g) and a human brown hair.The different variables that were used to implement the theoretical model for those two samples can be seen in Table 1.

Demonstration of the analytical model
The parameters for a thick silicone phantom containing TiO 2 and ICG (Table 1, first column) were used to visualize the model output at different steps of the simulation.The first behavior illustrated by the model is the change in OPL (Eq.5) which is a result of the cyclic variations in temperature caused by the photothermal laser and contrast agent.Our model simulates both the change in OPL at different depths inside the sample (Fig. 3a) and at a fixed depth over time (Fig. 3b).At this step, only the physical change in the sample is simulated.We still need to detect the change in OPL using an OCT system.This step is done using Eq.11, which calculates the experimental change in phase detected by a realistic OCT system.The experimental phase change is shown over time (Fig. 3c) at the same sample location depicted in Fig. 3b (z=400 μm).A Fourier transform is then performed on the experimental phase change signal (Eq.12, 13) and the PT-OCT signal is taken as the amplitude of the peak at the modulation frequency f 0 =500 Hz.Two characteristics of the signal can be observed from Fig. 3a) for ΔOPL and Eq. 6 for ΔΦ(z) : phase accumulation and shadowing.Phase accumulation describes the fact that ΔΦ increases with depth for a homogenous sample, and thus the photothermal signal will be detected as increasing, even though the amount of contrast agent in the sample is constant with depth.
Experimentally, this causes images to become brighter with depth even though the sample is homogeneous.Shadowing designates a region of high intensity in a region of the sample not containing any absorbers, simply because it is directly below a high absorbing region.Those two phenomena happen because, as seen in Eq. 5, when ΔT≥0, ΔOPL can only increase or remain stable with depth.ΔΦ is a function of optical path length which integrates with depth.Since we are never cooling the sample below its initial temperature, ΔOPL and ΔΦ can never decrease with depth, which causes phase accumulation and shadowing.Both phase accumulation and shadowing are detrimental to the image quality of PT-OCT.

Validation of the analytical model (monolayer)
A thick, one layer silicone phantom with 90 μg of ICG per gram of silicone was imaged to validate the result predicted by the theoretical model (parameters from Table 1, first column).
The power of the photothermal laser at the sample was varied between 0 mW and 16 mW.The resulting average PT-OCT signal can be seen in Fig. 4 for both the analytical model (blue) and the experimental validation (red).The two signals overlap throughout the field of view of the OCT system, for all powers of the PT laser.In each case, both the experimental and theoretical signal initially increases with depth, a result of phase accumulation.However, the signal becomes dominated by noise between 500 μm and 600 μm in depth, a result of the low reflectivity of the sample.This region in depth shows poor agreement between the theoretical signal and the experimental signal.However, this region of poor agreement corresponds to a region already dominated by noise (i.e.SNR< 0 dB).Thus, the model will not routinely be used to interpret data at these low SNR levels.Contrary to intuition, the noise floor (caused by the shot noise) is situated at 18 nm in PT-OCT signal (see Fig 4d).This is because the noise follows a 1/f behavior and happens to have a value of 18 nm at f=500 Hz.When the SNR approaches 0 dB the PT-OCT signal invariably increases or decreases to reach 18 nm.
Validation of the analytical model (multi-layer) and PT-OCT imaging depth A second experiment was used to investigate the maximum imaging depth of PT-OCT and validate the multi-layer model.The silicone phantom containing the human brown hair was imaged in such a way that the hair is at a different known depth (between 0 and 450 μm from the surface of the phantom to the center of the hair) every B-scan.The hair was considered to be 130 μm in thickness.An A-scan showing the experimental signal (Fig. 5a, red) from a hair at 80 μm in depth and the corresponding theoretical signal (Fig. 5a, blue) are shown to validate the multilayer analytical model.The model input to generate the analytical signal can be found in Table 1, second column.
To observe how the PT-OCT signal intensity decreases with the absorber depth, results from multiple hair depths were compared.For each hair depth, the PT noise was subtracted from the PT signal and this resulting signal was averaged over the cross section of the hair.The result was then recorded vs the hair depth (see Fig. We compared the evolution of the PT signal with depth and the evolution of the OCT signal intensity with depth (Fig. 5c).For this we measured the average OCT intensity of an area of scatterer inside the phantom, progressively choosing areas at greater depths inside the sample.
By comparing the point where the experimental PT signal and the experimental OCT signal have decreased by 50% of their initial value, it was determined that for this specific sample, the PT signal had decreased by half at 220 μm inside the sample, while the OCT intensity had decreased by half at 355 μm inside the sample.This indicates that PT-OCT has a smaller imaging depth than OCT, but of similar order.To test the robustness of the PT-CLEAN algorithm, a thin layer of phantom with TiO 2 and no ICG was added on top of the sample described previously and imaged with the same parameters (see Fig. 6c).The PT-CLEAN algorithm was applied with the same parameters and the resulting B-scan can be seen in Fig. 6d.Again, the PT-CLEAN signal is more constant over depth in the absorbing layer and close to the noise floor in the scattering layers.
An additional test for the PT-CLEAN algorithm is to differentiate between multiple concentrations of absorbers.Two thin layers of phantom were made with TiO 2 and respectively 11 μg and 45 μg of ICG per gram of silicone.Both layers were then mounted on top of a thick phantom with only TiO 2 with the high concentration layer on top (see Fig. 6e).The imaging parameters and the PT-CLEAN parameters were the same as for the previous two experiments.
The PT-CLEAN images can be seen in Fig. 6f.When the absorber concentration goes from high

Discussion
PT-OCT is an emerging technique for pre-clinical molecular imaging.It has been experimentally demonstrated in vitro, ex vivo and in vivo in the past and parts of its signal had been modeled.However, this is the first time that a comprehensive, customizable analytical model of a PT-OCT A-scan is presented for homogeneous and quasi-heterogeneous samples.We have validated our model against experimental data acquired in a homogeneous silicone phantom containing TiO 2 (scatterer) and ICG (absorber) at various photothermal laser powers, and in a multi-layer phantom containing TiO 2 and a human hair (absorber).Previous models of PT-OCT only described the general shape of the photothermal signal with depth and did not take into consideration the effect of noise.Our model is the first to match the data directly instead of being mathematically fitted to the experimental data points.Our model can also predict the A-scan obtained for different instrument design, imaging parameters, and sample composition.
By layering different sections of homogeneous samples, a heterogeneous sample was approximated.For in vivo drug delivery and molecular imaging experiments, it is expected that the contrast agent will diffuse through the sample and form a heterogeneous distribution.As a first step toward re-creating this situation, we created a three-layer phantom experiment and simulation: two large sections of scattering silicone surrounding a thin layer (130 μm) of absorber.A human brown hair was used as the thin absorption layer.This choice of absorber allowed us to easily vary the position of the absorbing layer in depth inside the sample and to assure that no leaking of the absorber into the scattering layer occurred.We measured the photothermal signal as a function of the hair depth and determined that the photothermal signal decreases with depth slightly more rapidly than OCT.This experiment also validated our model for a multi-layered sample, as the experimental values agreed with the theoretical values (NMSE=0.84).The maximum imaging depth for PT-OCT is strongly dependent on the type of sample being imaged, how many absorbing layers it has and where they are located.With our analytical model, we can now predict the imaging depth that will be obtained in practice for a given sample, which could not be done theoretically before.We can also select imaging parameters (e.g.photothermal laser power, photothermal modulation frequency, or the integration time of the detector) to maximize the signal strength in a specific sample.In the future, this model can be used to simulate a multi-layered sample that is a better approximation of a heterogeneous distribution of contrast agents inside a tissue sample.
Using our analytical model, we were able to simulate the PT-OCT signal from a homogeneous 7 μm thick absorber, which would be impractical to do without a model.This provided a point-spread function for our PT-OCT system.It was then possible to implement a modification of the CLEAN algorithm, which relies on the PSF of an instrument to improve image resolution and quality.The PSF should be generated with realistic optical parameters for the sample and instrument used to re-create the most accurate PT-OCT image.With the PT-CLEAN algorithm, we have demonstrated that phase accumulation and shadowing can be reduced in different layered samples and in a heterogeneous in vivo tumor sample, which is an improvement over traditional photothermal B-scans.Additionally, PT-CLEAN preserves the units of change in optical path length proper to PT-OCT and does not require any assumption on the signal shape or the sample composition, which allows its application to homogeneous and heterogeneous samples.In comparison with an algorithm proposed before [69], PT-CLEAN does not assume that the change in optical path length is proportional to the geometrical depth inside the medium (as seen in Fig. 3a, ΔOPL levels off with depth).It is therefore appropriate for complex, heterogeneous samples, and was thus the first algorithm of its kind to be applied to in vivo PT-OCT B-scans.Axial resolution for PT-OCT used to be difficult to define, since a small region of absorber could create a large shadow once the B-scan was processed.With the PT-CLEAN algorithm, the shadowing is eliminated and the effective axial resolution of PT-OCT is improved.
In conclusion, we have presented and validated the first analytical model of PT-OCT that can predict signal intensity with depth in homogeneous and multi-layered samples.We have also presented the PT-CLEAN algorithm to reduce phase accumulation and shadowing, effectively improving the axial resolution of PT-OCT. exp

Fig. 1 .
Fig.1.Flow diagram of the PT-CLEAN algorithm to remove phase accumulation and shadowing.A point spread function is created using the analytical model.The first pixel in depth with a PT-OCT signal above threshold is then selected to calculate the deconvolution kernel.Two new images are created; a B-scan from which the deconvolution kernel was subtracted and a new image where the selected pixel is recorded.The final image is reconstructed by convoluting each of the individual saved pixels to a 1-D Gaussian along the depth dimension.

Fig. 2
Fig. 2 PT-OCT instrumentation.The light from a superluminescent diode (SLD) is split between a reference and sample arm using a 50:50 fiber coupler (50:50).The photothermal excitation is provided by a Titanium:Sapphire laser (Ti:Sapphire) that is amplitude modulated by an acousto-optic modulator (AOM).HWP: Halfwave plate.PC: Polarization compensation.

Fig. 3
Fig. 3 Demonstration of the PT-OCT analytical output signal.(A) The change in optical path length produced in the sample by the photothermal effect is plotted as a function of depth and at a fixed depth z=400 μm over repeated Ascans (B).(C) The change in optical path length is detected as a noisy change in phase by the OCT system.The noiseless change in phase is overlaid in magenta to guide the eye.(D) A Fourier transform of the phase signal is taken to isolate the PT signal and PT noise.

Fig. 4
Fig. 4 Validation of the PT-OCT analytical model over an A-scan in a single layer homogeneous sample.The theoretical signal (blue) was generated from the model based on the optical properties listed in Table 1 (5mm ICG phantom).The experimental signal (red) is the average A-scan obtained experimentally by performing PT-OCT on a silicone phantom containing ICG as an absorber.The photothermal laser power at the sample was fixed at (A) 5.8 mW, (B) 9.2 mW and (C) 14.1 mW to show the behavior of the signal in relation to the noise floor.The noise floor can be seen in (D) when the photothermal laser is turned off.NMSE: Normalized Mean Square Error.

2 (
5b, red).The same value was generated by the analytical model (Fig. 5b, blue).The experimental and theoretical values overlapped for most of the hair positions.The normalized mean square error (NMSE = theoretical and experimental values at matched depth to quantify the agreement between the model predictions and the experimental result.(NMSE=0.84,perfect agreement when NMSE=1).

Fig. 5
Fig. 5 Validation of the PT-OCT analytical model over an A-scan in a multi-layered heterogeneous sample.(A) Theoretical (blue) and experimental (red) A-scan of a human brown hair 80 μm deep inside a scattering silicone phantom.(B) Average theoretical and experimental photothermal signal for different hair depths.Each hair depth was produced by a separate simulation or recorded on a separate B-scan.The signal has decreased by 50% when the hair is 220 μm in depth inside the phantom.(C) Experimental average OCT intensity for a section of scatterer at different depths inside the sample.The signal has decreased by 50% when the scatterers are 355 μm in depth inside the phantom.

Fig. 6
Fig. 6 Demonstration of the PT-CLEAN algorithm and comparison to traditional PT-OCT B-scans.(A) Thin ICG phantom on thick scattering (TiO 2 ) layer and (B) resulting B-scan after performing the PT-CLEAN algorithm.The horizontal dashed white lines indicate the limit of each layer.The solid white line is the average signal intensity vs. depth.(C) Thin ICG layer between a thin and a thick scattering layer and (D) after application of the PT-CLEAN algorithm.(E) Layer with a high ICG concentration, over a low ICG concentration layer, both over a scattering phantom and (F) after application of the PT-CLEAN algorithm.Scale bar: 100 μm.