1. Introduction

The brain requires a continuous supply of oxygen to provide energy and maintain proper function [1,2]. Pathological disruptions to the brain could impair oxygen delivery and brain metabolism [3]. As such, partial pressure of oxygen (pO2) is a critical indicator of cerebral physiology for characterizing alterations in brain function under normal and diseased conditions [4,5]. Two-photon phosphorescence lifetime microscopy (2P-PLIM) has been developed for non-invasive measurement of oxygen partial pressure (pO2) with micrometer-scale spatial resolution in brain tissue and vasculature in animal models. 2P-PLIM measures pO2 by quantifying the lifetime of an oxygen-sensitive phosphorescent dye. Previously, an oxygen probe PtP-C343 was created [6] and used for brain pO2 measurements in living rodents under awake and anesthetized conditions [7–9]. Recently, an improved oxygen-sensitive probe, Oxyphor 2P, was developed for deeper and faster brain pO2 measurement [10], and it has been applied on intravascular and tissue oxygen measurement of mice brain under baseline conditions and in response to functional activation [11,12]. The technique is based on oxygen-dependent quenching of phosphorescence, in which oxygen quenches radiative relaxation and shortens the phosphor’s lifetime [13]. The relationship between the oxygen partial pressure and the phosphorescence lifetime can be quantified via the Stern-Volmer equation [14]. Lifetime can be determined by calculating the decay rate of the measured phosphorescence emission curves. The standard approach to compute fluorescence or phosphorescence lifetime from time-domain measurements is through non-linear curve fitting to a mono- or multi-exponential decay function [8,11,12]. The method provides reliable calculation of lifetime and pO2. The accuracy of all lifetime computation techniques is heavily dependent on the number of photons collected during an experimental measurement. Additionally, however, the nonlinear fitting method is computationally intensive and subject to inaccuracies that depend proper selection of model parameters [15,16].

The “polar plot” or “phasor” analysis method has recently grown popular as a computationally-simple, model-free alternative to nonlinear curve-fitting for analyzing fluorescence lifetime imaging (FLIM) data [17,18]. The Fourier-based method was initially demonstrated by Jameson et al. in 1984 to analyze frequency-domain lifetime measurements of heterogeneous fluorophore mixtures [19]. It became more widely adopted ∼20 years later, after Digman et al. and Redford et al. applied it to time-domain FLIM measurements [17,20,21]. The technique entails computing the terms g(ω) and s(ω), the intensity-normalized cosine and sine transform of the time-resolved fluorescence measurement, at a specific angular frequency, ω. Generally, the phasor method makes no assumptions as to how the time-resolved measurements should be transformed. Consequently, applying this model-free computation to imaging data provides a simple and elegant technique to visually characterize variations in fluorescence lifetime by reversibly mapping pixels from a fluorescence image to points on a 2D phasor plot depicting the lifetimes of constituent fluorescent species. On a 2D phasor plot for a specific frequency ω, g(ω) and s(ω) constitute a cartesian coordinate pair for each time-resolved measurement. When applied to frequency-domain based measurements, phasor coordinates are equivalently represented in polar form by (M, φ), corresponding to the measurements’ modulation ratio and phase delay [21,22]. The phasor method has been applied extensively for time-correlated single photon counting (TCSPC) FLIM studies to assess FRET-based biosensors, resolve mixtures of spectrally-similar fluorophores, and characterize tissue autofluorescence and metabolism [17,23–25]. It has more recently been applied to analyze sample heterogeneity in other imaging methods such as quantitative MRI [26,27]. We previously applied phasor analysis to in vivo cortical FLIM measurements of endogenous, reduced nicotinamide adenine dinucleotide (NADH) to distinguish significant changes in mitochondrial activity [28,29].

Here, we apply the phasor method to analyze phosphorescence lifetime measurements lasting several hundred microseconds. Using phosphorescence lifetime computed by the phasor technique, we quantify cerebral pO2 in the awake mouse cortex. To avoid disadvantages associated with nonlinear curve-fitting and potentially reduce computation time, we sought to utilize the phasor method to analyze 2P-PLIM measurements of the novel oxygen-sensitizer, Oxyphor 2P, measured in cortical microvasculature. Rather than qualitative visual characterization of lifetime variation on the polar plot; however, we utilized the computed phasor coordinates to quantify phosphorescence lifetimes and corresponding pO2. We capitalized on the phasor method’s computational simplicity, and we hypothesized that the method may be less sensitive to factors that limit accuracy for curve-fitting. Because our established techniques for PLIM acquisition and analysis vary significantly from more common TCSPC FLIM methods, we implemented modifications to the phasor computation algorithm to achieve agreement with our conventional curve-fitting techniques. We validated the phasor-based lifetime and pO2 computation of experimental PLIM measurements by comparing curve-fitting results. To investigate factors that impact lifetime computation accuracy, we used simulation-based analyses to evaluate the effects of the longer PLIM excitation regime and different numerical methods for integration. We also compared both methods’ sensitivities to simulated parameters such as number of photons, varied levels of measurement noise, selected frequency, and theoretical lifetime. Choi et al. demonstrated pO2 quantification using phasor analysis using frequency-domain PLIM measurements in solution and cell preparations generated by temporal focusing wide field two-photon microscopy [30]. The results presented here advance upon previous efforts such as Choi et al’s., and, to our knowledge, signify the first time the phasor-based method has been applied to time-domain phosphorescence measurements with durations lasting hundreds of microseconds and for in vivo quantification of oxygen partial pressure, a vital metric for assessment of cerebral function. Our results indicate that the phasor method is a viable and robust method to computing cerebral pO2 from in vivo PLIM measurements. For our experimental measurements, phasor-based analysis produced pO2 values within ∼0.4–4% of pO2 computed by curve-fitting. Our simulation-based results also indicate that the phasor method may provide more robust pO2 measures in micro-environments exhibiting lower oxygen environments such as extravascular tissue or hypoxic tissue islets.

2. Methods

2.1 In vivo PLIM measurements

The phasor method was adapted and applied to experimental PLIM measurements to compute intravascular pO2. PLIM measurements were recorded from cortical microvessels of awake mice using a custom-designed multimodal microscope [31]. Experimental methods were analogous to those described previously [11,12]. Briefly, no less than 4 weeks prior to imaging, animals underwent surgery to implant a cranial window over the somatosensory cortex region and affix an aluminum post to the skull. Mice recovered for approximately 1-2 weeks after surgery and were then trained to remain stationary and calm when mounted in a custom cradle with their heads fixed. All procedures were approved by the Massachusetts General Hospital Subcommittee on Research and Animal Care and followed the Guide for the Care and Use of Laboratory Animals. Prior to imaging, ∼100 µl of Oxyphor 2P was injected at ∼6 mg/ml into the vasculature via retro-orbital injection. PLIM measurements were collected from intravascular locations up to 300 µm below the cortical surface. For each measurement, the dye was excited with a 10 µs-long train of 950 nm, ∼120 fs-long laser pulses from a tunable laser operating at 80 MHz (Insight DeepSee, Spectra Physics). Emitted phosphorescence photons (λ_em: 758 nm) were individually detected by a PMT connected to a discriminator circuit. Photons were acquired at 50 MHz sampling rate for a total 300 µs per measurement (10 µs laser excitation + ∼290 µs decay). At each intravascular location, measurements were repeated 2000 times to accumulate sufficient number of photons for further data processing. Recorded photon traces were temporally-binned into 2 µs intervals (t: 0–298 µs).

The experimental PLIM measurements were initially processed with our conventional curve-fitting methods to compute phosphorescence lifetimes and pO2 [8,10–12]. These curve-fitting results served as a reference and accuracy standard for our phasor-based analyses of the experimental measurements. For initial calibration studies, the developers of Oxyphor 2P rejected the first 5 μs of data after the falling edge of the excitation pulse to avoid confounding effects introduced by their system’s Instrument Response Function [10]. To maintain consistency with their calibration, only the decay portions (t = 16 -298 µs) of our measurements were fitted to a mono-exponential decay function with constant offset using nonlinear least squares functions in MATLAB. The fitted lifetime value was used to compute corresponding pO2 using an experimentally-derived calibration formula [10].

2.2 Modifications to the phasor method for PLIM

Using the curve-fitting results as a reference, we applied the phasor technique to compute phosphorescence lifetimes and quantify the associated pO2 values of our experimental Oxyphor 2P measurements. The phasor technique entails computing g(ω) and s(ω), of the time resolved fluorescence measurement f(t) at a specific angular frequency ω, respectively:

For time-domain FLIM measurements, ω is conventionally computed from the repetition period, T, of the pulsed excitation laser, $\omega = \frac{{2\pi }}{T}$, or an integral harmonic thereof. The computed values g(ω) and s(ω) therefore represent the intensity-normalized, real and imaginary parts of the measurement’s discrete Fourier transform. Using $g(\omega )$ and $s(\omega )$ as coordinates, each measurement is visualized on a 2-dimensional “phasor plot.” Phasor coordinates for mono-exponential decays lie along a “universal semicircle” centered at (g,s) = (0.5, 0) with radius = 0.5 [17,18,20,21].

To compute phosphorescence lifetime, τ, we utilized the a priori knowledge that Oxyphor 2P phosphorescence follows primarily mono-exponential kinetics. For a mono-exponential luminescence decay with lifetime, τ, $I(t )= {I_0}\textrm{exp}\left( { - \frac{t}{\tau }} \right)$, the corresponding phasor coordinates take the forms:

From g(ω) and s(ω), we tested the accuracy of computing τ using both the “modulation lifetime” (τ_m) expression and the “phase lifetime” (τ_ϕ) expression [30,32–34]. As detailed in the Supplement 1 (Table S1 and Fig. S1), τ_ϕ consistently agreed better with our curve-fitting analyses. Consequently, all lifetimes were computed using the τ_ϕ expression for a mono-exponential decay:

After computing τ, pO2 was calculated from the experimentally-derived calibration formula, similar to the curve-fitting approach [10].

To maintain consistency and achieve agreement between the two analysis techniques, we implemented modifications to the phasor method. Conventionally for FLIM data, the phasor integral is computed over the full time-resolved fluorescence measurement, and a calibration standard accounts for confounding influences of Instrument Response Function (IRF) [17]. For consistency with our curve-fitting method, our phasor computations were performed only over the isolated decay portion of the measurements (t = 16–298 µs). Additionally, to account for a background noise offset within our photon-counting measurements (which arises from background excitation, imperfect detection hardware, and stray ambient light), we subtracted the mean value of the final 14 µs (t = 284–298 µs) for each measurement. As we verified with simulations, our modifications for “trimming” the phasor integral boundaries (to include only the measurement’s decay portion) [35] and background removal effectively accounted for the adverse influences of the broader excitation pulse, potential presence of shorter lifetime components, and any prospective artifacts from the system’s Instrument Response Function. Additional simulation results are provided in the Supplement 1 (Figs. S2 and S3).

2.3 Phasor-based computation at different angular frequencies

For our phasor-based lifetime computations, the fundamental frequency, ω₁, is determined by the duration of the decay. Higher-order harmonics are designated as integral multiples thereof, ${\omega _n} = \frac{{2\pi \cdot n}}{{284\; \mu s}}$. To assess how lifetime computation accuracy varied with selected angular frequency, we computed the lifetimes and pO2 of our experimental measurements for n = 1-5.

As detailed in the results section, analysis of our experimental data indicated maximal agreement between curve fitting and phasor data when phasor computations were performed at the 3^rd harmonic, ${\omega _3}$. To investigate why the 3^rd harmonic consistently yielded better agreement, we sought to characterize the influence of harmonic frequency on the precision of lifetime computation for mono-exponential decays. Investigators have previously identified an “optimal” angular frequency that for a given lifetime τ, variations in Euclidian distance between phasor points are reportedly maximal when phasors are calculated at frequency ${\omega _{optimal}} = \; \frac{1}{\tau }$ [20,36]. We therefore assessed accuracy of phasor-based lifetime and pO2 calculations at the first 5 integral harmonics, and for our simulation-based investigations described below, we also tested the purported optimal frequency. The results were compared with lifetime and pO2 values computed from standard curve-fitting.

2.4 Comparison on experimental data

For each experimental PLIM measurement, we characterized agreement between phosphorescence lifetimes computed from both curve-fitting and phasor analysis. As the relationship between phosphorescence lifetime and pO2 is highly nonlinear, we also assessed agreement between the resultant pO2 values. Absolute differences for computed lifetimes and corresponding pO2s are reported in the results section. Residual errors were visualized between the experimental data and the mono-exponential curves generated with lifetimes from both methods. Comparisons on experimental measurements motivated us to perform more extensive simulation-based studies to assess differences in performance between both analysis techniques.

2.5 Simulation-based assessment of measurement parameters

We simulated phosphorescence decays in MATLAB (Mathworks, Framingham, MA) to assess how parameters such as angular frequency, integral boundaries, number of photons, and background noise affected the accuracy of lifetimes and pO2 computed with the phasor method and nonlinear curve-fitting.

2.5.1 Simulation with symbolic math

Using the Symbolic Math toolbox, we analyzed how isolating the decay portion of our measurements (t: 16 – 298µs) affects g(ω) and s(ω). We refer to the natural impulse response of phosphorescence emission as “Em_natural (t),” and modeled it as a mono-exponential decay function. For these simulations, the theoretical decay lifetime (τ) was held at 15 µs to represent an oxygen level pO2 = 69.2 mmHg, well within physiologically-normal vascular pO2 range. The 10 µs laser excitation gate from our 80 MHz laser was modeled as the sum of 800 Dirac delta functions, each separated by ${T_0} = $12.5 ns. We denote the complete phosphorescence measurement as “Em_full (t)” which was modeled by convolving the natural decay response with the excitation impulse train. The “Em_full (t)” profile is therefore represented as the summation of 800 exponential decays, each of which is time shifted by i*T₀:

$A$ is the initial photon count of each exponential decay. Experimental data indicates that the total number of photons in a single decay measurement is approximately 50 [12], such that the initial photon was determined by $A = \frac{{50}}{{\mathop \sum \nolimits_{t = 16}^{298} \textrm{exp}\left( { - \frac{t}{\tau }} \right)}}$. Importantly, the model for “Em_full(t)” is valid only if excitation duration and power is limited to avoid oversaturation of excited molecules within the focal volume.

The decay portion of simulated phosphorescence Em_full (t: 16–298 µs) identified as “Em_trimmed(t).” Phasor values of g(ω) and s(ω) were computed for “Em_full,” “Em_natural,” and “Em_trimmed” as the real and imaginary parts of each signal’s Fourier transform, divided by its respective time integral. Note that s(ω) represents the absolute value of the imaginary component of each function’s intensity-normalized Fourier transform.

2.5.2 Numerical simulations

In practice, the computational accuracies of phosphorescence lifetime and corresponding pO2 are limited by multiple experimental factors. To compare sensitivities of phasor analysis and curve-fitting methods to factors such as number of measured photons, relative amounts of background noise, and the nonlinear relationship between lifetime and pO2 values, we analyzed the computational accuracies for lifetime and pO2 over a broad range of numerically-simulated PLIM conditions and at various frequencies, ω, in MATLAB.

We simulated time-resolved distributions of detected photons by prescribing each photon’s arrival time [37]. Because our experimental measurements include contributions from both Oxyphor 2P emission and unwanted background noise, we modeled our PLIM measurements as the sum of detected photons from these two independent processes. Specifically, we initially assigned values for phosphorescence lifetime (τ), total number of photons (Nt) and the ratio between the phosphorescence “Signal” photons and “Background” noise photons (SBR). We divided the total number of measured photons (Nt) into either signal photons or background noise photons, and we modeled each photon’s arrival time as a randomly-determined number from either an exponential distribution (for emitted “Signal” photons) or a uniform distribution (for “Background” noise photons). A histogram of the arrival times was formed using 2 µs binning intervals to generate the time-resolved profiles. Note that only the decay part of the phosphorescence signal was simulated.

To determine the most representative values of Nt and SBR, we analyzed 201 experimental measurements. Nt was computed as the total number of experimentally measured photons from t = 16–298 µs, and τ was assigned to the lifetime value computed by conventional curve fitting. To find the representative SBR, we simulated time-resolved measurement profiles over a range of SBR values from ∼100 to 0.01 and checked for agreement between simulated phosphorescence and the experimental data using the reduced ${\chi ^2}$, defined as:

${E_i}$ and ${S_i}$ is the experimental data and the simulation data, $\sigma _i^2$ is the variance of the experimental data, N is the number of timepoints, and p is the number of floating fit parameters [34]. SBR of the experimental data was determined as value that gives the minimal reduced ${\chi ^2}$. Based on the average values of total photon number and SBR across our experimental measurements, 96000 and 28 were used as representative values for numerical simulation of phosphorescence decay.

Using the numerical simulation of phosphorescence emission, we independently examined the influence of Nt, SBR, and τ for lifetime and pO2 calculations. As detailed in the results section, we determined the influence of varying each parameter independently while maintaining the other 2 parameters at fixed values.

3. Results and discussion

3.1 Analysis of experimental PLIM measurements

Using experimental PLIM measurements as test data, we modified the conventional phasor computation method for phosphorescence lifetime quantitation. As discussed in section 2.2, the modifications included (a) confining the phasor integral calculation to only the measurements’ decay portion (t = 16–298 µs) and (b) removal of undesired background offset from our photon counting measurements. After removal of the average background signal, negative photon counts were set to zero. Phasor-based computations were performed at ${\omega _n} = \; \frac{{2\pi \cdot n}}{{284\; \mu s}}$ for n = 1-5.

Preliminary characterizations revealed that the standard numerical methods for computing Eq. (1) and Eq. (2) for TCSPC FLIM data yielded insufficient accuracy for our PLIM measurements. With the 284 μs PLIM decays divided into 2 μs binning intervals, the relative width of integral partitions is almost twice as large compared to our FLIM measurements. As a result, the simple rectangular method for numerically computing definite integrals proved inadequate for PLIM data [38]. The inaccuracies were readily apparent when visualized on a 2D polar plot, as the phasor coordinates computed for theoretical mono-exponential decays from 2–38 µs did not align well with the ideal universal semicircle (Supplement 1). Integration accuracy was restored after implementing improved quadrature methods (Simpson’s rule).

Figure 1(a-c) display representative phosphorescence lifetime data computed from an experimental PLIM measurement using both conventional curve-fitting and phasor analysis, respectively. Figure 1 b shows the phasor representation of phosphorescence lifetime computed using the fundamental harmonic. Phasor plot generated using higher harmonics are shown in Fig. S5. As seen in Fig. 1.b, phasor points computed from Oxyphor 2P phosphorescence lie close to the unit semicircle, indicating that phosphorescence kinetics follow a mono-exponential decay. A zoomed-in phasor plot can be found in Supplement 1 to illustrate Oxyphor 2P’s mono-exponential kinetics. We analyzed residuals between the experimental data and theoretical mono-exponential decays. The decays were generated using computed lifetimes from both curve-fitting and phasor analysis at the fundamental frequency, ω₁. Example data recordings, theoretical reconstructions, and corresponding residuals are provided in Figs. 1 c-e. In general, the residuals for both methods appear randomly scattered around zero at all time points. The residual plots indicate that the analyzed segment of Oxyphor 2P phosphorescence can be modeled well by a mono-exponential function, and furthermore, the phasor method and curve-fitting method perform comparably for computing the mono-exponential PLIM lifetime. We also tested bi-exponential model for fitting the measurement data and found no significant improvement. Examples of bi-exponential fit and the associated residuals are provided in Supplement 1.

 

Fig. 1. (a) Two-photon intensity image of field of view (FOV). Colored circles denote selected intravascular locations for phosphorescence lifetime measurements and are color-coded by the corresponding pO2 values. (b) Polar plot representation of phosphorescence lifetime of the selected pixels. Filled circles represent pO2-coded phasors of experimentally measured Oxyphor 2P phosphorescence. Hollow circles represent phasors of theoretical mono-exponential decays, all of which fall along the unit semicircle (c) Mono-exponential fits of measured phosphorescence photon emissions with lifetimes given by non-linear curve fitting and phasor analysis. (d-e) Weighted residuals of the exponential fits. (f) Difference between lifetime computed by phasor analysis and curve fitting method of all data sets: $\mathrm{\Delta }\tau = {\tau _{phasor}} - {\tau _{curve{\; }fitting}}$. First five harmonics (${\omega _1}\sim {\omega _5}$) were used in phasor approach for calculating lifetime. (g) Difference between pO2 computed by phasor analysis and curve fitting method of all data sets: $\mathrm{\Delta }pO2 = pO{2_{phasor}} - pO{2_{curve\; fitting}}$.

Download Full Size | PDF

Figures 1(f) and (g) display the mean and standard deviation of the difference between curve-fitting and phasor analysis for computed lifetimes and pO2, respectively, across all 201 analyzed experimental measurements. The differences are provided for phasor-derived values computed at ${\omega _n}\; ({\; n = 1,2,3,4,5} )$. In general, the differences of lifetimes computed with phasors at the first five harmonics were within 1 µs of the lifetimes computed with curve-fitting, while the corresponding phasor-derived pO2 values showed larger variations. When analyzing our experimental data, we consistently found that τ and pO2 computed at the 3^rd harmonic agreed best with the curve-fitting based results, yielding mean pO2 difference of 0.01 ± 1.22 mmHg (0.38 ± 2.60%), while the 1st harmonic yields pO2 difference of -2.22 ± 2.21 mmHg (-4.44 ± 3.24%) relative to curve-fitting results.

Conventionally for frequency-domain measurements, emission lifetimes are computed using both the phase and modulation of the measured emission signal relative to excitation. For time domain measurements, phasor coordinates can be utilized to compute both the “phase lifetime” and the “modulation lifetime” (Supplement 1 Eq. (1) and Eq. (2)) [30,32–35]. Table S1 shows the lifetime values calculated using the phase and modulation (${\tau _\phi }$, ${\tau _m}$) of the measured decay for pixels in Fig. 1 a. Fig. S1 shows the phase and modulation lifetimes and their associated pO2 values for all 201 Oxyphor 2P measurements. Our experimental results indicate “phase lifetime” and corresponding pO2 agree much better with results obtained from conventional nonlinear curve-fitting. “Modulation lifetimes” were consistently overestimated, which, thereby underestimated pO2.

3.2 Symbolic math simulations

We used simulations to verify that our modifications to the phasor analysis algorithm account for the PLIM method’s longer excitation paradigm, the IRF, and the potential influence of shorter-lifetime components. The Symbolic Math toolbox in MATLAB permitted analysis over a broad range of ω while avoiding discretization errors associated with conventional numerical methods. Figure 2 a illustrates the initial 50 μs of the temporal profiles, where “Em_full” represents phosphorescence emission during and after excitation with radiative lifetime of 15 µs, “Em_trimmed” represents “Em_full(t : 16–298 µs)”, and “Em_natural” is a mono-exponential decay function that represents the impulse response of Oxyphor 2P emission. Figure 2 b and 2 c display g and s for each profile as a continuous function of ω. Notably, g(ω) of “Em_full” diminishes more sharply with increasing ω, while s(ω) of “Em_full” first increases, crests, and then diminishes more sharply. In practice, a 10 µs excitation pulse train is necessary to generate sufficient detectable phosphorescence signal for our experimental measurements. As Fig. 2 illustrates, however, the longer excitation regime substantially alters both the temporal profile and frequency content of measured Oxyphor 2P emission. The profiles of “Em_trimmed” overlap completely with those of “Em_natural”, verifying that our method of limiting the phasor integral to the decay portion effectively removes the influence of longer excitation gate.

 

Fig. 2. Symbolic simulation of phosphorescence lifetime decay. (a) truncated temporal profiles of simulated theoretical decay, PLIM measurement, and trimmed portion of PLIM measurement. (b-c) Corresponding phasor coordinates g and s, as a continuous function of the angular frequency $\omega $, for each profile vertical lines show the position of the first 10 harmonics. $1/\tau $ is reportedly the “optimal frequency” for maximal phasor separation and is displayed as the black vertical line.

Download Full Size | PDF

Figure 2(b) and (c) illustrate how the amplitudes of phasor coordinates g(ω) and s(ω) vary with frequency for a mono-exponential decay function. The vertical lines indicate corresponding integral frequency harmonics at which phasor coordinates are often calculated and visualized for time-resolved measurements. The bold vertical line corresponds to the purported “optimal frequency” for lifetime calculation [20,36]. Local maxima for g(ω) and s(ω) occur at ω = 0 and near the third harmonic ω₃, respectively. Lifetime, τ, is proportional to the ratio of the two coordinates, and theoretically, its value should be constant across all frequencies. The profiles suggest how the sensitivity of lifetime calculation varies across ω in practice. At higher harmonics where both g(ω) and s(ω) are smaller, discrepancies in computed g and s values are more probable due to Poisson variability and background noise and would yield significantly larger errors in computed lifetimes. Consequently, quantification of lifetimes or other related metrics using higher harmonics (n > 4) is inadvisable due to reduced accuracy.

Additional simulation results are provided in the Supplement 1. Using methods developed by Laine et al. [35], Figs. S2 and S3 illustrate the effects of “trimming” the integral boundaries on lifetime computation accuracy over a broad range of lifetimes. Similar to Laine et al’s. findings with shorter lifetimes, our results show that lifetime computation accuracy is greatly diminished when the lifetimes exceed a threshold percentage of the integration period (∼13%). Our procedure for trimming the integration period does not adversely affect lifetime computation accuracy within the calibrated lifetime range of Oxyphor 2P. However, accuracy is diminished for lifetime values above or below the calibrated range. Symbolic simulation results using longer and shorter lifetimes are included in Fig. S8. Relative to τ = 15 µs, the profiles for g(ω) and s(ω) at τ = ∼12 µs and τ = ∼31.2 µs appear slightly broadened and constricted, respectively. The optimal frequency lies closer to ω₄ for τ ∼ 12 µs (pO2 ∼ 100 mmHg) and closer to ω₁ for τ ∼ 31.2 µs (pO2 ∼ 10 mmHg).

3.3 Numerical simulations: influence of background noise

In practice, the presence of background noise distorts lifetime calculation, both for curve-fitting and frequency-domain calculations. Radiative photon emission follows Poisson-like statistics, and therefore, produces substantial variability for photon counting measurements. Furthermore, artifacts from imperfect detection electronics, leakage from the excitation laser, and detection of ambient light introduce randomly-distributed noise that further confounds photon-counting measurements. The instrument background noise introduces an offset to the recorded photon profiles, and we typically remove its influence by subtracting the mean value of the final 10 -15 µs from PLIM measurements before fitting. To determine the influence of background noise, we performed numerical simulations to assess the errors caused by varying signal-to-background ratio on lifetime and pO2 calculation accuracy for curve-fitting and phasor analysis.

Time resolved PLIM measurements were modeled by simulating the detection times for a discrete number of photons. The relative proportion of “Signal” photons and “Background” photons varied by adjusting the Signal to Background Ratio (SBR) from 1–100, while the total number of photons (Nt) remain fixed at 96000.

Simulations were repeated 50 times for each value of SBR. The relationship between measured lifetime τ and corresponding pO2 is highly nonlinear. Therefore, to assess whether the influence of SBR varies at different oxygenation levels, we performed these simulations with 2 different theoretical lifetime values τ, representative of measurements from highly oxygenated arterioles (pO2 = 100 mmHg, τ = 12.0 µs) and from venules with lower oxygen levels (pO2 = 30 mmHg, τ = 22.7 µs).

Figure 3 displays how computed lifetime and pO2 varied from their “ground-truth” values at different values of SBR. For each value of SBR, PLIM measurements were simulated, and the lifetimes and pO2 values were computed by curve-fitting and by phasor analysis at the first 5 integral harmonics and $\omega = \frac{1}{\tau }$. Mean differences in lifetime from ground truth are generally below 1µs for all SBR values under simulated high and low oxygen settings, and notably smaller lifetime differences are observed in the high oxygen, short lifetime setting. As Fig. 3 c and d illustrate, the small error in lifetime corresponds to discrepancies in pO2 calculations of 1-3 mmHg. Variability in pO2 was found to be larger at higher levels of oxygen, as expected from the Stern-Volmer relationship. Interestingly, increasing levels of SBR from 1 to 100 yielded only modest improvements to the computation accuracy for both lifetime and pO2. The results suggest that our practice for background removal effectively minimizes potential calculation errors caused by instrumentation noise over a wide range of experimental conditions. Our simulations suggest no appreciable differences in accuracy for curve-fitting and phasor computation at harmonics n=1-5.

 

Fig. 3. Computed lifetime and pO2 accuracy from simulated PLIM data at different signal to background ratios (SBRs). Simulations were performed for high (arteriole simulation, pO2 = 100 mmHg) and moderately low (venule simulation, pO2 = 30 mmHg) intravascular oxygenation levels.

Download Full Size | PDF

3.4 Numerical simulations: influence of number of photons

Using numerical simulations, we assessed how the number of measured photons affects lifetime and pO2 calculation accuracy for the phasor method and curve-fitting. Time resolved PLIM measurements were modeled as described in sections 2.5.2 and 3.3. For these analyses, SBR was held fixed at 28, representative of our experimental data, while the total number of photons (Nt) was varied from 1000–100000. Simulations were repeated 50 times for each value of Nt.

Differences between ground truth lifetime and pO2 and their computed values are displayed in Fig. 4 for different collected photon totals. Again, to account for the nonlinear relationship between τ and pO2, we performed these simulations with 2 different theoretical lifetime values τ, representing highly oxygenated arterioles (pO2 = 100 mmHg, τ = 12.0 µs) and venules with lower oxygen levels (pO2 = 30 mmHg, τ = 22.7 µs). As expected, both phasor analysis and curve-fitting show pronounced sensitivity to the number of collected photons. As indicated by the amplitude of the error bars, the uncertainty in computed pO2 values ranged from 4 to more than 20 mmHg when the total number of collected photons was less than 10000. Increasing the number of photons above 25000 yielded notable improvements to the measurement accuracy for both lifetime and pO2. As the radiative photon emission process follows Poisson statistics, increasing the total photon counts reduces the statistical uncertainty. The variability introduced by this Poisson-like radiative process impacts our calculation accuracies much more substantially than SBR [16]. Our results indicate that approximately 50000 - 90000 photons are necessary to achieve pO2 accuracies within ±2 mmHg for our experimental measurements for pO2 quantitation by both curve-fitting and phasor analysis. For our experiments, this corresponds to ∼1000–2000 repeated 300 µs measurements at each location [12]. The simulated results also suggest that none of the tested frequencies had superior performance for phasor calculations, contrary to the analyses of our experimental data. Instead, we observed generally comparable accuracies between curve-fitting and phasor-based results for all tested frequencies.

 

Fig. 4. Computed lifetime and pO2 accuracy from simulated PLIM data for different numbers of total collected photons. Simulations were performed for high (100 mmHg) and moderately low (30 mmHg) intravascular oxygenation levels.

Download Full Size | PDF

3.5 Numerical simulations: variability across the physiological range of pO2

Lifetime and pO2 follow a highly nonlinear, Stern-Volmer relationship for phosphorescent oxygen sensors, where high pO2 levels correspond to shorter lifetimes, and vice versa [14]. As a result, even for a constant level of variability in computed lifetime accuracy, pO2 accuracy still varies disproportionately across different oxygenation levels. Using numerical simulations, we analyzed how the lifetime and pO2 calculation accuracy change over the physiological range of pO2 for both phasor method and curve-fitting. Time resolved PLIM measurements were modeled as described in section 3.3. For these analyses, SBR was maintained fixed at 4 and Nt was fixed at 50000 to represent a measurement in deep scattering tissue with the minimum acceptable photon total and maximum acceptable background noise. Simulations were performed for over the lifetime range of 9.5 to 35 µs, corresponding to pO2 range of 140–5 mmHg. 50 simulated photon traces were tested for each value of τ.

Figure 5 shows differences (mean ± std dev) between theoretical values and our simulated PLIM measurements across the physiologically-relevant range of pO2. Computations of τ and pO2 were performed using curve-fitting and phasor method at the first harmonic and at $= \; \frac{1}{\tau }$ . The shaded regions represent standard deviation, and their magnitudes indicate that variability in lifetime measurement was reasonably low (< 1 µs) over the entire pO2 range, with lower variability at higher pO2. As seen in Fig. 5 b, despite lower variability in computed τ, the variability in pO2 increased with higher oxygen levels up to ±2 mmHg. The results indicate reasonable levels of uncertainty for τ and pO2 under these simulated, sub-optimal conditions. In practice, however, fewer photons would be emitted under shorter-lifetime, higher oxygen conditions, and the variability would be predictably higher in arterioles than veins. Again, both curve-fitting and phasor analysis yielded comparable levels of computation accuracy. Interestingly, Fig. 5 b indicates substantially lower uncertainty for phasor analysis compared to curve-fitting at low oxygenation levels (pO2 < 40 mmHg), while curve-fitting appears to perform better at higher oxygenation levels (pO2 > 100 mmHg). Additional investigation is required; however, these results suggest that phasor analysis could yield better accuracy for pO2 measurements in tissue or under hypoxic conditions, while curve-fitting could be more robust for measurements under hyperoxic conditions.

 

Fig. 5. Computed lifetime and pO2 accuracy from simulated PLIM data with theoretical pO2 levels spanning normal tissue physiology. Differences between ground truth and computed values are displayed for curve-fitting and phasor analysis at the first harmonic and $\omega = 1/\tau $.

Download Full Size | PDF

4. Conclusion

This study indicates that the Fourier-based phasor analysis method shows promise for quantifying partial pressure of oxygen from PLIM measurements. Fourier-based calculations have previously been applied for characterizing heterogeneous distributions of phosphorescence lifetimes recorded with frequency-domain measurements [39]. Additionally, pO2 was previously computed with phasor analysis of frequency-domain PLIM measurements of a ruthenium-based oxygen sensor in solutions and cell preparations generated by temporal focusing wide field two-photon microscopy [30]. The oxygen-sensitive ruthenium sensor emits phosphorescence with lifetimes ranging from 400–600 ns. To our knowledge, the present study is the first to apply the phasor technique to time-domain phosphorescence measurements lasting hundreds of microseconds, and the first compute oxygen partial pressure in vivo. This report presents the first phasor-based observations of spatially-resolved pO2 measurements performed in vivo, collected from the microvasculature of the awake mouse cortex. Experimental- and simulation-based analyses show that the phasor method generally yields comparable performance to conventional curve-fitting methods for quantifying phosphorescence lifetime and pO2. For both analytical methods, numerous parameters influence computation accuracy, and among the factors that we investigated, the total number of recorded photons had the largest influence on computation accuracy. Future investigations will focus on the feasibility of phasor analysis to improve measurement accuracy and computation speed in extravascular brain tissue and under low-oxygen conditions.