Optical property recovery with spatially-resolved diffuse reflectance at short source-detector separations using a compact fiber-optic probe

: We describe a compact fiber-optic probe (2 mm outside diameter) that utilizes spatially-resolved diffuse reflectance for tissue optical property recovery. Validation was performed in phantoms containing Intralipid 20% as scatterer, and methylene blue (MB), MnTPPS, and/or India ink as absorbers. Over a range of conditions, the reduced scattering coefficient was recovered with a root mean square error (RMSE) of 0.86-2.7 cm − 1 (average error = 3.8%). MB concentration was recovered with RMSE = 0.26-0.52 µM (average error = 15.0%), which did not vary with inclusion of MnTPPS (p = 0.65). This system will be utilized to determine optical properties in human abscesses, in order to generate treatment plans for photodynamic therapy.


Introduction
Spatially-resolved diffuse reflectance has been widely applied for recovery of optical properties from biological tissue [1][2][3][4]. This method allows for separation of the effects of tissue absorption, represented by the absorption coefficient (µ a ), and scattering, represented by the scattering coefficient (µ s ). In highly scattering media, the reduced scattering coefficient (µ s '), which includes the effects of scattering anisotropy, is often used to characterize tissue scattering. The goal of these prior studies has typically been to determine optical properties that can distinguish between diseased and healthy tissue [5,6], or allow for treatment planning and monitoring of light-based therapies such as photodynamic therapy (PDT) [7,8] or photothermal therapy [9]. These treatment planning solutions rely on forward models of light propagation that can be used to optimize placement of light sources and treatment parameters to maximize efficacy and minimize risk [10,11].
The current study is envisioned in the context of treatment planning for PDT of deep-tissue abscess cavities. Abscesses develop in response to bacterial and/or fungal infection, and result in the formation of a fibrous pseudocapsule enclosing a purulent fluid collection. While percutaneous drainage and antibiotic prescription have become standard of care, cure rates remain low in some cases [12] and recurrence or complication following drainage remain problematic [13,14]. For these reasons, we have initiated a Phase 1 clinical trial examining the use of methylene blue mediated PDT at the time of percutaneous abscess drainage (ClinicalTrials.gov Identifier: NCT02240498). In a prior simulation study, we demonstrated that the portion of the abscess wall achieving a desired fluence rate threshold decreased with increasing absorption at the wall of the abscess [15]. This, along with other studies showing the importance of accurate optical property values in PDT treatment planning for hollow cavities [16], motivates the need for measurements of optical properties at the abscess wall. The clinical scenario also constrains the maximum size of the diffuse reflectance probe, as the entire assembly is required to be inserted through the 12 French (F) catheter typically used for abscess drainage. This limits the outer diameter of the probe to approximately 2 mm, with source-detector separations constrained by fiber construction and the thickness of the sheath surrounding the optical fibers.
Many surface-contact, spatially-resolved diffuse reflectance probes have been created for recovery of tissue optical properties. A number of these utilized approximations to the radiative transport equation in order to perform this extraction, with a large amount of work focused on the diffusion approximation [1,17,18]. In cases where the assumptions of the diffusion approximation cannot be met, other investigators have explored the application of the P 3 approximation to the radiative transport equation [19][20][21]. While this approximation improves accuracy at shorter source-detector separations and is applicable for a wider range of optical property combinations, there are still fundamental limitations on probe geometry and allowable combinations of optical properties. One particular concern is validity at short source-detector separations, which may be required in cases where the size of the fiber-optic probe is restricted by clinical considerations.
For this reason, many modern solutions have moved to Monte Carlo (MC) based optical property recovery [22][23][24][25]. In this case, a MC model of the optical probe is created and a lookup table is generated by simulating detected diffuse reflectance at a range of optical property combinations. Absorption and scattering can then be determined by comparison to the lookup table, while constraining possible values either spectrally or spatially. For example, Hennessy et al reported accurate optical property recovery using a MC lookup table and assumptions of the absorbers present in the sample [22].
However, these MC solutions are not without limitation. As described above, many implementations require knowledge of the expected absorption spectra present in the sample [22,26]. This is sufficient for many cases, where oxy-and deoxy-hemoglobin are known to be the main absorbing species. For our clinical application, optical properties have never been quantified in human abscesses. While it is reasonable to expect that hemoglobin will be a major contributor to absorption, multiple bacterial species have been shown to produce porphyrins with absorption in the visible spectrum [27,28], a fact which has been exploited to perform red-light inactivation of these strains [29]. Therefore, we require the ability to extract optical properties in the presence of unknown absorbers. For this reason, we focus on spatially-resolved data at individual wavelengths in a compact fiber-optic probe, with additional exploration of the added benefit of knowledge of underlying absorption spectra.
Here we describe the design and validation of a surface-contact, spatially-resolved diffuse reflectance instrument that will allow for recovery of tissue optical properties from the interior of human abscess cavities. This incorporates a compact fiber-optic probe with an outside diameter (OD) of 2 mm, and a Monte Carlo (MC) based optical property recovery algorithm. Results are provided for a wide range of tissue-simulating phantoms, including cases with multiple absorbing species.

Fiber optic probe and spectroscopy system
The optical probe used for diffuse reflectance measurements consists of nine individuallyaddressable optical fibers (see Fig. 1(a)), each with a 200 µm core diameter and 0.22 numerical aperture (NA). One fiber acts as a source for spectroscopy, while the other eight function as detectors. These fibers are encapsulated in a flexible assembly with an outer diameter of 2 mm. This outer diameter is sufficiently small to pass through the lumen of the 12F drainage catheters (Flexima APDL, Boston Scientific, Marlborough, MA) typically used for the percutaneous drainage procedure described in the Introduction. At the distal end of the probe, each of the fiber faces is arranged in a "plus" pattern, as shown in Fig. 1(a). With a designed spacing of approximately 300 µm between adjacent fiber centers, this results in actual source-detector separations ranging from 288 µm to 1.30 mm. The probe was manufactured by Pioneer Optics Company (Bloomfield, CT), and the entire assembly is designed to withstand high-level disinfection with Cidex OPA (Advanced Sterilization Products, Irvine, CA). Upon receipt of the probe, the detector fiber positions relative to the source fiber were imaged with a digital camera, as shown in Fig. 1(b). Images were captured with a ruler in frame, in order to calibrate distance. Eight detector fibers were selected for the current probe design based on a balance between adequate sampling of radial distance, practicality of device manufacture, and robustness of the final product. Spectroscopy measurements were performed with this probe using the optical system design shown in Fig. 1(c). For reflectance measurements, the source is a fiber-coupled broadband tungsten halogen lamp (HL-2000-HP-FHSA, Ocean Optics, Inc., Largo, FL). For fluorescence measurements, the source is a fiber-coupled laser diode at 640 nm (OBIS FP 1193841, Coherent, Inc., Santa Clara, CA). Note that while the system is capable of performing fluorescence measurements, only the reflectance functionality is reported here. Based on the desired measurement, the proper source is selected with a 2 × 1 multimode optical switch (FSM 1 × 2, Piezosystem Jena, Inc., Hopedale, MA) and routed to the source fiber in the optical probe. Probe detection fibers are connected to an 8 × 1 multimode optical switch (FSM 1 × 8, Piezosystem Jena, Inc.) for individual read-out. The output of the detection switch is connected to a free-space collimator (RC02SMA-P01, Thorlabs, Inc., Newton, NJ), which routes the detected signal through a motorized filter flip mount (MFF101, Thorlabs), before being re-focused into the fiber with an identical collimator. This flip mount allows for a fluorescence emission filter (LP02-664RU-25, Semrock, Rochester, NY) to be rotated into the beam path for fluorescence measurements. The final output is detected by a high-sensitivity, thermoelectrically-cooled spectrometer (QE Pro, Ocean Optics). The system is controlled by an encrypted, password-protected laptop computer using a custom interface built in LabVIEW (National Instruments, Austin, TX).

Data collection and correction
For each of the experimental conditions described in section 2.5, reflectance spectra were captured sequentially at each of the eight detector fibers using the optical switches described above. In order to utilize the dynamic range of the spectrometer, the integration time was set automatically for each detector fiber at the time of measurement. This was done by making an initial measurement with an integration time of 50 ms, and dividing the peak signal by a value corresponding to 66% of the maximum value allowed by the spectrometer. This multiplier was then used to determine the final integration time for each fiber, which ranged from 40 ms to 5000 ms.
Each spectral measurement was corrected for dark background, integration time, lamp power, and system throughput. For background correction, dark measurements were taken separately from light measurements for each fiber and source using the same integration time, room lighting, and probe orientation but without the system lamp enabled. This signal was subtracted from the measured signal. Background subtracted spectra were subsequently divided by integration time in ms and measured lamp power in mW, with lamp power measured using a calibrated integrating sphere power measurement system (Labsphere, North Sutton, New Hampshire). To evaluate the system response, the probe was placed with its distal end flush with the input port on an integrating sphere (Labsphere, North Sutton, New Hampshire) and a NIST traceable lamp (LS-1-CAL, Ocean Optics, Largo, FL) was used to illuminate the sphere. The measured spectra were background subtracted, divided by the known lamp spectrum, and normalized to a maximum value of one to generate individual fiber system responses. The background-subtracted spectra captured for experiments were divided by these system response spectra to correct for sources of spectral distortion within the system. Prior to taking each set of experimental data, a calibration measurement was taken to account for possible changes in fiber throughput and fluctuations in source power. To characterize this response, the probe was situated with its distal end flush with the input port on the integrating sphere, and measurements were taken using the probe's source fiber and system lamp. This calibration was performed for each detector fiber sequentially, corrected using the procedure described above, and smoothed using the Savitzky-Golay method across 21 data points with a second-degree polynomial. Experimental data from phantom measurements were then divided by the corresponding calibration data for each detection fiber. As described in the Introduction, the clinical application for this spectroscopy system requires that the probe maintain high level disinfection prior to use. Therefore, the probe could not be simply mounted above a flat diffuse reflectance sample with a fixed mount, as others have done [30,31], as contact with a non-sterile mount would compromise the probe disinfection. An integrating sphere was used in place of a diffuse reflectance standard due to reduced sensitivity to probe height and angle compared to a flat standard (data not shown), and the ability to perform non-contact calibration.

Monte Carlo model
In order to recover optical properties from spatially-resolved diffuse reflectance measurements, a Monte Carlo (MC) model of the fiber optic probe was created. This was created as an addition to the publicly available CUDAMCML [32], as we are interested in running graphics processing unit (GPU) accelerated simulations for a semi-infinite medium. The probe was modeled as a 2 mm diameter cylinder with a refractive index of 1.46, with individual fibers having a radius of 100 µm, refractive index of 1.46, and NA=0. 22. The source and detector fibers were positioned on the probe face identically to the measured positions described in section 2.1. Photon packets were launched from the source fiber within the NA, and allowed to propagate within tissue as normal. At the interface between the probe and tissue, photon packets were first checked for transmission or reflection based on calculation of Fresnel coefficients using the refractive indices of the probe and tissue. Photon packets that were transmitted into the probe within the radius and NA of a given detector fiber were scored as detected with their remaining photon weight. Photon packets that were transmitted into the probe but did not strike a detector fiber within the NA were terminated. For photon packets striking the air-tissue interface outside of the diameter of the probe, Fresnel coefficients were calculated using the refractive indices of air and tissue.
Simulations were performed with the probe model positioned vertically with its distal face in contact with the tissue boundary in order to replicate experimental conditions. A lookup table was generated by running MC simulations at all combinations of µ a = 0.0001-10 cm −1 and µ s = 1-250 cm −1 , with n=1.37 and using the Henyey-Greenstein phase function with g=0.7. Scattering anisotropy was chosen based on the formula reported by van Staveren et al at 665 nm [33], while refractive index was based on values reported by Dong et al [34]. Absorption coefficients were sampled in increments of 0.0001 cm −1 from 0.0001-0.001 cm −1 , 0.001 cm −1 from 0.001-0.01 cm −1 , 0.005 cm −1 from 0.01-2 cm −1 , 0.1 cm −1 from 2-5 cm −1 , and 0.5 cm −1 from 5-10 cm −1 , while scattering coefficients were sampled uniformly in increments of 0.25 cm −1 . Detection values at each detector fiber were combined across optical property combinations to generate a 3D MC library for comparison to measured data. Surface plots of this MC library are shown in Fig. 2 for detectors at source-detector separations of 288 and 974 µm. All simulations were run with 10 8 photon packets to reduce noise, compared to fewer photon packets, on a Quadro RTX6000 GPU with 24 GB of GPU memory (NVIDIA Corporation, Santa Clara, CA). For 10 8 photon packets, this results in variance of ∼1% between repeated identical simulations, as compared to ∼16% for 10 6 photon packets.

Recovery of optical properties
Measured, corrected data were used to recover optical properties by comparison to the 3D MC library described above. This was done for two cases: (1) fitting of data at each wavelength independently, and (2) simultaneous multi-spectral fitting using the known absorber and scatter spectral shapes. For both cases, it was first necessary to calibrate measured data to the simulation library, as the scales are different (i.e. counts/mW/ms vs. simulated photon weight). In order to do this, fully corrected data for methylene blue phantoms at n=36 combinations of absorber and scatterer concentration (phantom condition 2, described below) were divided by simulated data with identical optical properties at each detector fiber to generate ratio spectra, where D n (λ,d) is the detected spectrum for the nth phantom condition at detector d, Sim n (λ,d) is the simulated spectrum at the same optical properties as the nth phantom condition at detector d, and R n (λ,d) is the ratio spectrum for the nth phantom condition at detector d. For each detector fiber, these ratio spectra were then averaged across phantoms and wavelength to generate a single ratio value for each detector fiber, where R n (λ,d) are the ratio spectra described in Eq. (1) for the nth phantom condition and detector d, n is the number of phantom conditions, n λ is the number of wavelengths, and R avg (d) is the ratio value for detector d. This averaging was performed to reduce the sensitivity of the calibration to variation in one specific phantom. Data for all phantoms were divided by the ratio values at the appropriate detector fibers to generate a calibrated, detected signal, before comparison to the MC library.
For fitting at individual wavelengths, the fitting algorithm looped over individual wavelengths and calculated the difference between the measured data and the MC library for each detector fiber independently. The error metric was then calculated at each combination of (µ a ,µ s ') as, where MC(µ a ,µ s ',d) is a matrix of simulated photon weights at all combinations of (µ a ,µ s ') for detector d, and D cal (d) is the detected signal, divided by the MC calibration ratio described above, at detector d. The recovered optical properties at a given wavelength are then determined by finding the combination that minimizes the error function. This process is repeated to build up recovered spectra, µ a (λ) and µ s '(λ), with no explicit association between values at adjacent wavelengths. In order to convert the recovered absorption spectra to concentrations of known chromophores, the error function, was minimized using constrained, non-linear optimization, where µ a,recovered (λ) is the spectrum of absorption coefficients recovered at each wavelength, conc(A) is the concentration of absorber A, and B(A,λ) is the absorption basis function for absorber A.
For the case of multispectral fitting, absorption was assumed to be a superposition of known absorption spectral shapes and scattering was assumed to be of the form, where λ 0 is a normalization wavelength set to 665 nm. Using initial guesses for absorber concentrations, µ s '(λ 0 ), and b, the error function, was then minimized using constrained, non-linear optimization, where Sim(conc,µ s '(λ 0 ),b,λ,d) is the simulated detected spectrum for detector fiber d at the specified combination of absorber concentrations, conc, and scattering parameters and D(λ,d) is the measured spectrum at detector fiber d. In this way, measured spectra at all detector fibers are used simultaneously to determine absorption and scattering spectra. Interpolation was used to determine Sim(conc,µ s '(λ 0 ),b,λ,d) at combinations of (µ a ,µ s ') that were not explicitly simulated. For both cases, fitting was performed over the wavelength range of 500-740 nm, as this is the range over which the chosen chromophores have appreciable absorption. Data were re-sampled at a wavelength resolution of approximately 1 nm, in order to improve fitting speed. All optical property recovery was performed in MATLAB (R2019b, The MathWorks, Inc., Natick, MA), with constrained non-linear minimization performed using the fmincon function and the active-set algorithm. In all cases, absorption and scattering were constrained to be positive, and the value of b in Eq. (6) was constrained to be in the range of 0.5-4. For each set of spectra, extraction of optical properties took approximately 3 seconds for individual wavelength fits, and approximately 1 second for full spectral fitting on an Intel i7-9700 CPU with 16 GB of memory.

Phantom preparation and experimental conditions
Three categories of phantoms were used throughout these experiments. The first was a combination of Intralipid 20% as scatterer (Fresenius Kabi AG, Bad Homburg, Germany) and India Ink as absorber (Higgins No. 4418, Chartpak Incorporated, Leeds, MA), used to determine the system's ability to recover absorption or scattering at a single wavelength. Second was a combination of Intralipid 20% as scatterer and Methylene Blue as absorber (MB, American Regent, Inc, Shirley, NY), used to determine the system's ability to detect full absorption and scattering spectra of known shape. The third consisted of a combination of Intralipid 20%, with both MB and Mn(III) meso-Tetra (4-sulfonatophenyl) porphine (MnTPPS, frontier scientific, Logan, UT) as absorbers, which was used to determine the system's ability to extract individual absorption spectra in the presence of background absorption. MB was selected as an absorber due to its usage as a photosensitizer in the clinical application described in the Introduction. MnTPPS was used due to the similarity in spectral shape to that of hemoglobin, ease of preparation, and stability. MnTPPS has successfully been implemented in previous phantoms in our lab [35,36].
The first set of phantoms consisted of 1L water, Intralipid as scatterer, and India ink as absorber. Large volumes were used to replicate a semi-infinite phantom. Intralipid was added to produce separate phantoms with µ s = 25, 50, 75, and 100 cm −1 at 665 nm. For each of these, India ink was added to generate µ a =0.05, 0.1, 0.2, 0.5, 0.75, and 1.0 cm −1 at 665 nm. This resulted in 24 unique phantoms.
The second set of phantoms consisted of 1L water, Intralipid 20% as a scatterer, and MB as an absorber. Intralipid volumes were added to produce scattering coefficients of µ s = 25, 50, 75, 100, 125 and 150 cm −1 at 665nm. To each volume of Intralipid increasing volumes of MB were added to produce absorption coefficients of µ a = 0.05, 0.1, 0.2, 0.5, 0.75 and 1.0 cm −1 at 665nm, corresponding to MB concentrations of 0.32µM to 6.4µM. This resulted in 36 unique phantoms consisting of all combinations of these absorption and scattering coefficients.
The third set of phantoms consisted of 1L water, Intralipid 20%, MB, and MnTPPS. The volume of Intralipid was determined to achieve a scattering coefficient of 100 cm −1 . This was held constant and the volume of MnTPPS stock was added to obtain absorption coefficients of µ a = 0.1, 0.2 and 0.5 cm −1 at 562nm from MnTPPS, corresponding to concentrations of 4.5µM to 22.5µM. For each constant volume of Intralipid and MnTPPS, MB was added to achieve absorption coefficients µ a = 0.05, 0.1, 0.2, 0.5, 0.75 and 1.0 cm −1 at 665nm from MB alone. This resulted in 18 unique phantoms that allowed for analysis with the presence of multiple absorbers.
To determine the necessary volumes of Intralipid for each phantom, Intralipid 20% was diluted with distilled water to a range of 0.02-0.37%. The spectra were collected using a commercial absorption spectrophotometer (Cary 50 Varian, Santa Clara, CA) with a 1 cm path length cuvette. Assuming negligible absorption, the measured optical density was considered reflective of single scattering effects and treated as its scattering spectra, as we have done previously [37], resulting in spectra comparable to those reported by Michels et al [38].These spectra were used to fit a power-law curve of the form OD(λ)= aλ −b . The spectra were normalized to 665 nm and fit simultaneously to the power law-curve in order to find b. Next, the relationship between concentration and magnitude of the spectra was used to find the scaling factor a. Due to potential variance between bags of Intralipid, the same batch of Intralipid was used wherever possible and new scattering spectra were taken for each unique stock.
India ink absorption was quantified using the method described above over concentrations of 0.6-4.9 µL/mL water. We found that the absorption spectra of India ink was fit well (R 2 =0.99) with an expression of the form µ a (λ) = c · [ink] · λ −d , where c=1963, d=1.23, and [ink] corresponds to ink concentration in units of µL ink per mL solvent. To determine the necessary volume of MB, solutions ranging from 1.5-11.6 µM were added to a 1cm path length cuvette and absorption spectra were taken using the same spectrometer described above. To find the shape of the basis spectrum, the data were normalized to 665 nm, the average optical density was taken across each wavelength, and the data were smoothed and divided by their maximum value. The scaling factor was found by taking the relationship between concentration and the magnitude of individual spectra. This same method was used for MnTPPS over a range of 4.8-36.2µM and normalizing the basis to the maximum value at 562nm.
To create the phantoms, a 3L container was spray painted black (Rust-Oleum Flat Black, Vernon Hills, IL) and filled with 1L of distilled water. Intralipid was added and another set of spectra were taken. Next MnTPPS was added if it was used and another set of spectra were taken. Lastly, India ink or MB was added in increasing concentrations. The phantom was constantly stirred using a magnetic stir bar, also painted black. Once MB was introduced, the phantom preparation was performed in the dark to avoid potential bleaching effects under room lights. The probe was held in surface contact with the phantom using laboratory clamps while measurements were taken.

Statistical analysis
Recovered optical properties and chromophore concentrations are represented throughout by mean ± standard deviation across phantoms with the same relevant optical property (e.g. across all Intralipid concentrations with [MB] = 6.4 µM). Recovered absorber concentrations using the two recovery methods described in section 2.4 were compared to known values and each other using repeated measures one-way ANOVA. Pairwise comparisons were performed using Tukey's test. Bland-Altman analysis was used to examine bias between the two optical property recovery methods, with the resultant bias compared to a theoretical mean of 0 using a one-sample t-test. Differences between recovered and known phantom values are also reported as percent error compared to known values and root mean square error (RMSE). All statistical analyses were performed using GraphPad PRISM (v6, GraphPad Software, Inc., San Diego, CA).

Ink phantoms
Data collected in phantoms containing India ink were fit using the procedure described in section 2.4 at 665 nm. Results are shown in Fig. 3, with recovered µ a and µ s ' values averaged across phantoms at identical ink or Intralipid concentrations, respectively. µ s ' was not found to be significantly different between known and recovered values (p=0.082), with RMSE = 2.7 cm −1 and average error of 6.5%. While recovered µ a increases monotonically with ink concentration, recovered absorption was systematically higher than the known (Bland-Altman bias = 0.085 cm −1 (-0.23-0.40 cm −1 95% confidence interval (CI), p=0.011), with RMSE = 0.15 cm −1 and average error of 23.6%. This error was particularly pronounced at the lowest ink concentration examined.

Methylene blue phantoms
Representative data after full correction as described in section 2.2 are shown in Fig. 4. Here, color corresponds to detector fiber, while line style corresponds to the MB concentration. Radial fall-off and the impact of increasing MB concentration are readily apparent.
These data were fit using the MC method described above, both at each wavelength individually and over the full spectra simultaneously using the known shape of MB absorption. Recovery of µ a spectra for both of these cases are shown in Fig. 5, with summaries provided in Fig. 6 for both absorption and scattering.
Both methods showed monotonic increases in recovered MB concentration with known MB concentrations, with RMSE=0.52 µM for individual wavelength fits (average error = 7.0%) vs. RMSE=0.47 µM for the full spectral fits (average error = 20.6%). However, MB concentrations returned by the full spectral method were significantly higher than the known concentration (p=0.0024), whereas the other method returned values similar to known values (p=0.75). As with the ink phantoms, relative error was worst at the lowest MB concentration examined and relative error decreased with increasing MB concentration. The reduced scattering coefficient at 665 nm was recovered very well, with RMSE = 0.95 cm −1 for individual wavelength fits (average error = 3.1%) and RMSE = 0.86 cm −1 for full spectral fits (average error = 2.8%). For both cases, µ s ' was not found to be significantly different from known values (independent: p=0.75; full spectrum: p=0.46).    6. Recovery of (a) MB concentration and (b) µ s ' at 665 nm using either independent wavelength or full spectral fits. Symbols represent mean values across phantoms with identical optical properties for each method, with error bars corresponding to standard deviation. Solid red lines indicate perfect agreement. Note that error bars are included for all cases, but are not always visible.

Methylene blue and MnTPPS phantoms
Representative corrected data are displayed in Fig. 7. The effects of MnTPPS absorption can be seen around its Soret band and in the region from 500-625 nm. Here, phantoms consisted of constant Intralipid and MnTPPS concentration, and the MB concentration was escalated. These data were again fit using the previously described optical property recovery algorithms. MB concentration was fit well with both methods, as shown in Fig. 8(a), with RMSE = 0.26 µM for independent wavelength fits (average error = 18.2%) and RMSE = 0.32 µM for full spectral fits (average error = 14.4%). Recovered concentrations were similar to known concentrations in both cases (independent: p=0.14; full spectrum: p=0.94), with values also similar between fitting methods (p=0.096). When compared to phantoms that only included MB, the recovered MB concentration was very similar in the presence of MnTPPS (p=0.65). Data collected in the presence of MnTPPS tended to return a slightly lower MB concentration (Bland-Altman bias = -0.16 µM (-1.8 µM -1.5 µM), though this difference was not significant (p=0.65). This indicates that concentrations of relevant chromophores are recovered accurately, even in the presence of relatively large background absorption (up to 0.5 cm −1 ). MnTPPS concentration was also recovered well, as shown in Fig. 8(b). Error was slightly higher than for MB recovery, with RMSE = 1.36 µM for individual wavelengths (average error = 21.5%) and RMSE = 1.44 µM for spectral fits (average error = 14.9%). Recovered [MnTPPS] was significantly higher than known values for independent wavelength fits (p=0.003), though values recovered by the two recovery algorithms were not significantly different from each other (p=0.19).

Discussion
Here, we demonstrate a diffuse reflectance spectroscopy system incorporating a compact fiberoptic probe that can be inserted through the 12F catheters used for routine abscess drainage. Using a Monte Carlo based recovery algorithm, we showed accurate recovery of both absorption and reduced scattering coefficients at a single wavelength, as well as for full detected spectra. Critically, the accuracy of these extracted optical properties was not affected by the presence of multiple absorbers. For full spectra data, optical property recovery was performed both with and without knowledge of the absorbers present in the sample. We found that recovered chromophore concentrations were similar between the two methods, indicating that optical property recovery in the presence of unknown absorbers is feasible and accurate. This is vital for our clinical application, as this will represent the first measurement of tissue optical properties in human abscess cavities.
As described in the Introduction, many other research groups have demonstrated extraction of tissue optical properties from spatially-resolved diffuse reflectance measurements. The main advantages here are the utilization of short source-detector separations, reduced restriction of viable optical property combinations, and lack of reliance on knowledge of underlying absorbers present in the sample. While approaches utilizing the diffusion approximation are widespread, this approximation relies on relatively large source-detector separations and high transport albedo [a'=µ s '/(µ a +µ s ')]. Typically, diffusion-based techniques are limited to source-detector separations greater than approximately 1 mean free path, and transport albedos greater than 0.98 [39]. For these reasons, the P 3 approximation has also been explored [19]. While this approximation improves the range of validity, there are still fundamental limits on minimum source-detector separation, with a typical minimum being on the order of 500 µm [19]. These restrictions are not present in the current study, where we utilize source-detector separations as short as 288 µm and transport albedos as low as 0.89. However, it should be noted that systems utilizing analytical approximations for inversion of spatially-resolved diffuse reflectance measurements can produce more precise determinations of absorption, due to the larger optical paths through tissue detected by these systems. While the reduced scattering coefficients recovered here are highly accurate and reproducible, recovery of absorption is relatively worse, due to the short path length through the sample.
The accuracy of our results is comparable to those reported by other groups using MC lookup table methods. For example, LaRiviere et al demonstrated errors in recovery of µ a and µ s ' ranging from 8.4-65.9% and 6.4-19.9%, respectively [23]. In phantom validation, Voulgarelis et al showed an average error of 1.9% in recovery of µ a and 4.4% for µ s ', although this was only performed across a relatively high µ a range of ∼4-9 cm −1 [24]. Naglič et al demonstrated average errors of 15.6% in absorption and 5.8% in reduced scattering at source-detector separation ranging from 220-1100 µm [25]. A major limitation of this study is that recovery was only accurate for µ s ' values greater than 15 cm −1 . Hennessy et al showed excellent recovery of both hemoglobin concentration and µ s ' (2.4% and 1.7%, respectively), although their experimental setup required a large spectrograph and camera [22]. By comparison, we show an average error of 15.0% in µ a and 3.8% in µ s ', calculated across all measured phantom conditions. We demonstrate similar or better recovery of µ s ' compared to these previous studies, but comparable or worse recovery of µ a . However, these other methods often require physically larger, inflexible probes [6,25,31], which are not suitable for our clinical application, or assumptions on the absorbers present in the sample [22]. As our clinical application will represent the first measurements of optical properties in human abscess cavities, for which there are not prior data on the absorbers expected to be present, constraining the absorbers during data inversion would fundamentally limit the knowledge gained from these measurements.
As described in the Introduction, our primary clinical application will be in measuring tissue optical properties at the wall of human deep-tissue abscess cavities. The described system will be incorporated into our ongoing Phase 1 clinical trial of methylene blue PDT at the time of abscess drainage (ClinicalTrials.gov Identifier: NCT02240498). In this trial, methylene blue is infused into the abscess at a concentration of 1 mg/mL, allowed to incubate for 10 minutes, and then aspirated before administration of treatment light. As MB is highly absorbing, it will be important to measure abscess wall optical properties both before and after MB incubation, in order to quantify native abscess absorption and scattering, as well as MB uptake during the incubation period. These factors have been studied extensively in relation to oncological applications of PDT [40], and we will apply this same level of rigor to antimicrobial PDT. While photosensitizer concentration is a critical component for efficacious PDT, along with light dose, a threshold MB concentration has not yet been determined for antimicrobial applications. Indeed, a recent review reported a wide range of MB concentrations (0.3-60 mM) that have proven efficacious in clinical antimicrobial PDT studies [41]. In in vitro studies, cells are typically incubated at a particular MB concentration and rinsed before illumination, without quantification of MB uptake [42][43][44]. For in vivo studies, MB is generally applied topically, again without quantification of MB uptake [45][46][47]. As a result, there is a need to link the actual concentration of MB taken up by cells to response to PDT, with this work ongoing. Further, for topical application, MB uptake may be depth dependent, as has been shown for other photosensitizers [48][49][50][51]. This would affect the assumption of homogeneous optical properties that we currently employ, and would have an impact on the deposited light dose. As we have access to multiple source-detector separations, we may be able to perform coarse examination of this depth dependence, as has been demonstrated in multiple studies [52][53][54].
For these reasons, our current focus is on the effect of measured optical properties on the delivered light dose. In the clinical setting, we will use the recovered optical properties to generate patient-specific PDT treatment plans to optimize light dose, as we have previously demonstrated for head and neck cancer applications [55]. The reliance of delivered light dose (fluence rate) at the abscess wall on absorption has been previously demonstrated [15], with ongoing work demonstrating large increases in patient eligibility for abscess PDT based on generation of individualized treatment plans (manuscript in preparation). In the present study, we demonstrate an average error of 15% in the determination of µ a . Let us consider the case of an abscess with an actual value of µ a =0.5 cm −1 at the wall. Based on our results in preparation, this will require 311 mW to be delivered to the cavity to achieve a fluence rate of 4 mW/cm 2 at 95% of the abscess wall for a representative abscess. If we impose the current 15% error on this µ a value, resulting in an assumed absorption coefficient of 0.58 cm −1 , the average optical power to achieve the light dose target increases to 335 mW. A 15% error in µ a therefore results in a 7.7% error in the delivered optical power. If our goal is to achieve an error in the prescribed optical power of less 10%, we can therefore tolerate approximately 20% error in the extraction of µ a , which has been achieved here.
While the present system has demonstrated accurate results and will be utilized for our ongoing clinical study, there is still room for improvement. Particularly for the envisioned multi-center Phase 2 trials, simplicity and ease of use will be paramount. As described briefly in the Methods section, the present system also incorporates the capability for fluorescence measurements. While this potentially allows for recovery of photosensitizer concentrations at a lower range than that achievable by reflectance, it does add a significant amount of complexity to the system design. This includes a laser source for excitation, an optical switch for source selection, and a motorized filter for fluorescence emission. Future studies will focus on the utility of fluorescence measurement, in order to determine whether this capability is necessary for additional prototypes.
The MC based recovery algorithm also depends on precise measurements of source-detector separation at the probe face and generation of a large MC lookup table for each detector fiber. While the MC code utilized is GPU accelerated and executed on a high-performance GPU, this still represents a significant time expense. Further, this MC library would need to be generated for each probe employed. For this reason, among others, many groups have shifted to machine learning based optical property recovery [56][57][58]. Rather than using analytical approximations or large MC libraries, these techniques are fully data-driven. This can result in accurate optical property recovery for relatively simplistic designs, without the need for cumbersome model generation. For example, Nguyen et al demonstrated that a variety of machine learning models outperformed MC-based methods in both accuracy and run-time using simulated diffuse reflectance spectra [59]. We plan to explore these methods in the context of our application, in order to make optical property recovery simpler and more robust.
We acknowledge some limitations in the present study. While we endeavored to include a wide range of optical properties, extreme values of µ a and µ s ' were not included. Therefore, the absolute range of validity is not known, though transport albedo ranged from 0.89-1 in the generated phantoms. Optical properties were also assumed to be homogeneous, both radially and axially from the probe axis. As described above, methylene blue may have depth dependent uptake, which would result in violation of this assumption. Future design iterations will focus on multi-layered MC models incorporating depth-dependent optical properties, as discussed above. Additionally, we acknowledge the relative complexity of the present spectroscopy system compared to previously reported systems. As described above, future iterations on the design will stress simplicity and ease of use, compared to the present prototype. Finally, the MC-based optical property recovery described here is dependent upon an accurate model of the probe, which requires precise measurement of source-detector separations on the probe face, as well as calibration between the physical system and the MC model. Incorporation of alternative probe designs would require similar quantification of relevant quantities, as well as generation of new MC models, look-up tables, and calibration ratios between the MC model and experimental values. As mentioned above, machine learning approaches may alleviate this time-consuming process, and will be explored in greater detail. Further, periodic re-calibration with reference phantoms should be performed to account for any temporal change in the optical system and probe.