1 Introduction

Molecular oxygen (O2) has a central role in cellular respiration where it is used in the production of adenosine triphosphate (ATP); the main molecule supplying energy to metabolic processes for supporting cell and organ function [15]. The primary site of oxygen delivery from blood to organ cells is the microcirculation. Therefore, microcirculatory oxygenation is considered a parameter of key (patho)physiological importance [26, 28]. Hence, monitoring microcirculatory oxygenation under (patho)physiological conditions has been the focus of many studies, for example, to assess the severity of organ hypoxia during shock and to guide treatment strategies aimed at improving organ oxygenation [7]. Methods for quantitative assessment of microcirculatory oxygen tension (pO2) include invasive techniques, such as oxygen electrodes [3, 32], and non-invasive techniques, such as phosphorimetry [2, 21, 24, 29, 30, 35].

Phosphorimetry is a non-invasive optical method for measuring oxygen tension quantitatively in vivo, based on oxygen-dependent quenching of phosphorescence of metallo-porphyrins, such as palladium porphyrin [4, 29, 30]. The phosphorescence decay trace is measured and fitted using a mono-exponential curve. The pO2 is calculated using the Stern–Volmer relationship, which directly relates the phosphorescence decay rate (or lifetime) [25]. Under (patho)physiological conditions, however, it has been well established that microvascular oxygenation could be spatially heterogeneously rather than homogenously distributed within the phosphorimetric measurement volume and that, within one organ, hypoxic tissue could co-exist with normoxic tissue [7, 36]. These hypoxic tissue fractions could potentially lead to organ dysfunction and therefore their detection is essential. Instead of mono-exponential curve fitting of the phosphorescence decay trace, providing a single pO2 value for the phosphorimetric measurement volume, multi-exponential curve fitting could be used for the analysis of the decay traces, which provides a distribution (or histogram) of pO2 values [1113].

The fundamental method for multi-exponential curve fitting analysis of fluorescent and phosphorescent decay traces for recovering pO2 distributions is the exponential series method (ESM) [9]. An alternative analytical approach, the maximum entropy method (MEM) [16], was compared to the ESM and both methods were shown to be equally powerful in resolving underlying distributions, but since the ESM was significantly faster than the more complex MEM, ESM has been suggested as the method of choice [23]. Furthermore, the MEM was recently questioned as recovered pO2 distributions sometimes contained nonphysiologically high pO2 values. This was later attributed to the reduced accuracy of the MEM for high pO2 values due to noise [27, 33]. Vinogradov and Wilson have reported a deconvolution method that used similar algebraic methodology as in the ESM, but where function minimization was approached though the application of quadratic programming principles [31]. This method provided results similar to the ESM in a significantly reduced computation time. Refinement of the ESM was proposed by Golub et al. by transforming the measured decay traces using a heterogeneity factor, as introduced by Zheng et al., which accounts for the difference between decay traces obtained from continuous quencher concentration distributions (i.e., the measured decay trace) and the equivalent discrete distributions (i.e., the decay trace obtained from the pO2 histogram as described by the ESM) [6, 34].

Although it is generally accepted that the oxygen-quenched phosphorescence and delayed fluorescence decay traces can be analyzed using multi-exponential curve fitting, its application until now has been limited to a few (patho)physiological studies, probably because the reliability of the recovered pO2 histograms has never been extensively evaluated and lacks documentation. The aim of the present study is to extensively evaluate the ability of multi-exponential fitting analysis to adequately determine pO2 histograms from simulated decay traces at different signal-to-noise ratios (SNRs). In the first part of the study, we tested whether the application of the heterogeneity factor as proposed by Golub et al. actually improves the recovery of pO2 histograms as this has never been shown; how noise affects the recovery of pO2 histograms; and whether the effects of noise could be suppressed by application of a median filter. However, the adverse effect of the heterogeneity factor is an amplification of both the signal and the noise at the end of the decay trace. Therefore, rather than transforming the phosphorescence decay trace using the heterogeneity factor and thereby amplifying the noise in the latter part of the trace, here, we included the heterogeneity factor in the fit model as a modified version of the method developed by Golub and Zheng et al. and applied by Mik and Johannes et al. [6, 19, 20, 34]. In the second part of this study, we applied multi-exponential curve fitting analysis for the recovery of bimodal pO2 distributions, with a hypoxic and a normoxic mode, at several SNRs and tested the ability of the analysis procedure to detect hypoxia. In the third part of the study, we applied multi-exponential curve fitting analysis for the recovery of microvascular pO2 histograms measured in the rat kidney in a model of endotoxemic shock and fluid resuscitation.

2 Materials and methods

2.1 Generation of decay traces

Analysis software for the simulation of pO2 distributions, conversion to decay traces, multi-exponential curve fitting analysis, and recovery pO2 histograms was written in LabView (version 8.6, National Instruments, Austin, TX, USA).

The continuous input pO2 [mmHg] distributions simulated for this study were described by a single or double Gaussian function f(pO2), where f(pO2) ≥ 0 and ∑f(pO2) = 1. The simulated pO2 distributions were subsequently converted to continuous decay rate [k (μs−1)] distributions via the Stern–Volmer equation: k(pO2) = k 0 + k q  · pO2 [25]. Here, k 0 (μs−1) and k q (μs−1 mmHg−1) are the phosphor-dependent decay rate at pO2 = 0 mmHg and quenching constant, respectively. Decay traces P(t) were obtained by substitution of the pO2 distributions into:

$$ P\left( t \right) = \int {f\left( {{\text{pO}}_{2} } \right)e^{{ - \left( {k_{0} + k_{q} \cdot {\text{pO}}_{2} } \right) \cdot t}} {\text{dpO}}_{2} } $$
(1)

To the decay traces, having an amplitude of 1 at t = 0 as described by the Gaussian function, variable uniform white noise (representing mainly thermal noise) was added to set the SNR (defined as the maximum phosphorescence amplitude/peak-to-peak noise) from 100 to 10, representing SNRs typically found in experiments (a typical decay trace is presented in Figs. 1b, 5a).

Fig. 1
figure 1

a The heterogeneity factor as described by Golub et al. [5] for transformation of the measured decay trace. b A typically measured decay trace (black), having a signal-to-noise ratio (SNR) of ~20, and the corresponding transformed decay trace (gray)

2.2 Multi-exponential curve fitting

The generated decay traces, corresponding to continuous pO2 distributions, were analyzed by multi-exponential curve fitting using a constrained Levenberg–Marquardt procedure, a procedure which minimizes the weighted sum of the squared residuals (i.e., data-fit) χ2. Multi-exponential fitting was performed by standard fitting of a sum of exponentials to the decay trace (i.e., the ESM) [9].

In the ESM, the decay trace is fitted by a sum of N exponentials with fixed, uniformly spaced decay rates, k i (pO2), and variable pre-exponential fraction coefficients, f i (pO2), which are to be recovered. The histogram bins have a width of 2δ mmHg and are centered on pO2 = δ + 2iδ mmHg, where i = 0, 1,…N − 1:

$$ P\left( t \right) = \sum\limits_{i = 0}^{N - 1} {f_{i} \left( {\delta + 2i\delta } \right)e^{{ - \left[ {k_{0} + k_{q} \cdot \left( {\delta + 2i\delta } \right)} \right] \cdot t}} } $$
(2)

Equation 2 does not take into account the difference between decay traces obtained from continuous quencher concentration distributions (i.e., the measured decay trace Eq. 1) and the equivalent discrete distributions (i.e., the decay trace obtained from the pO2 histogram as described by the ESM). To correct for this issue, Golub and Zheng et al. proposed to multiply the decay traces with a time-dependent heterogeneity factor [6, 34]:

$$ P\left( t \right)^* = P\left( t \right)\left[ {{\frac{{e^{{k_{0} t}} k_{q} \delta t}}{{\sinh \left( {k_{q} \delta t} \right)}}}} \right] = \sum\limits_{i = 0}^{N - 1} {f_{i} \left( {\delta + 2i\delta } \right)e^{{ - k_{q} \cdot \left( {\delta + 2i\delta } \right) \cdot t}} } $$
(3)

However, the heterogeneity factor amplifies both the signal and the noise at the end of the decay trace, and therefore, in this study, rather than transforming the phosphorescence decay trace using the heterogeneity factor we included the heterogeneity factor in the fit model and thereby prevented amplification of the noise in the latter part of the trace (Fig. 1). The fit model is then described by:

$$ P\left( t \right) = \left[ {{\frac{{\sinh \left( {k_{q} \delta t} \right)}}{{e^{{k_{0} t}} k_{q} \delta t}}}} \right]\sum\limits_{i = 0}^{N - 1} {f_{i} \left( {\delta + 2i\delta } \right)e^{{ - k_{q} \cdot \left( {\delta + 2i\delta } \right) \cdot t}} } $$
(4)

2.3 Simulations

In this study, the decay trace according to Eq. 1 is simulated with a pO2 range from 0 to 300 mmHg with a resolution of 1 mmHg. Simulated sampling rate was set to 1 MHz (1 sample/μs). k 0 and k q were 1/251 and 273 μs−1 mmHg−1, respectively. Measurement time was set to 2/k 0, i.e., 502 μs. These values correspond to values previously found for the applied phosphorimeter and phosphorescent dye as described below.

Tolerance for χ2 for the fit optimization was 10−3. The initial fit coefficients were set at f i  = 1/N, where N is the number of bins in the pO2 histogram, so that the sum of all fraction coefficients equals 1. As a constraint, the fraction coefficients were not allowed to be scaled below 0. The histogram bin width was set to 10 mmHg (δ = 5) and the histogram pO2 range was set from 0 to 150 mmHg (N = 15).

In the first set of simulations, we tested whether the application of a noise-reducing median filter or the heterogeneity factor to the ESM would actually improve the recovery of pO2 histograms from (noisy) decay traces, having an SNR of 100 and 10. In the second set of simulations we applied the ESM for the recovery of bimodal pO2 distributions, with a hypoxic and a normoxic mode at different ratios, at several SNRs and tested the ability of the analysis procedure to detect hypoxia.

2.4 Data analysis

Data plotting was performed in GraphPad Prism (GraphPad Software, La Jolla, CA, USA). To allow clear visual comparison of pO2 distributions recovered using the different analytical approaches, the mean + SD values in the pO2 distributions (n = 50 runs per simulation) were plotted and connected with lines. In the first set of simulations, the recovered pO2 histograms were compared to the input distributions using the Kolmogorov–Smirnov test, similarly to Golub et al. [6]. The D value provided by this test was used to lookup a corresponding P value (also depending on the number of histogram bins), where P > 0.05 indicates that the input and recovered histograms are not significantly different (i.e., are similar). In the second set of simulations, the fraction of the first bin of the recovered pO2 histograms was plotted versus the simulated degree of hypoxia and analyzed using linear regression analysis to evaluate the adequacy of the detection of hypoxia in bimodal pO2 distributions.

2.5 In vivo experiment

To demonstrate the (patho)physiological relevance of multi-exponential curve fitting analysis of in vivo obtained phosphorescence decay traces for recovering microvascular pO2 histograms, measurements were performed in the rat kidney (n = 5) at baseline, during endotoxemic shock, and after fluid resuscitation. The protocol for this experiment has been approved by the animal ethics committee of the Academic Medical Center at the University of Amsterdam and has been described elsewhere [11].

Rats were anesthetized with an intraperitoneal injection of 90 mg/kg ketamine (Nimatek®; Eurovet, Bladel, the Netherlands), 0.25 mg/kg medetomidine (Domitor; Pfizer, New York, NY, USA) and 0.05 mg/kg atropine-sulfate (Centrafarm, Etten-Leur, the Netherlands). After tracheotomy the animals were mechanically ventilated with a FiO2 of 0.4. Body temperature was maintained at 37.0 ± 0.5°C by external warming. The ventilator settings were adjusted to maintain an arterial pCO2 between 35 and 40 mmHg. A polyethylene catheter (outer diameter = 0.9 mm; Braun, Melsungen, Germany) was placed in the right carotid artery to monitor arterial blood pressure and heart rate.

The left kidney was decapsulated and immobilized in a Lucite kidney cup (K. Effenberger, Pfaffingen, Germany) via a ~4 cm incision of the left flank. The phosphorescent dye (Oxyphor G2, λexcitation = ~440 nm, λemission = ~800 nm; Oxygen Enterprises, Philadelphia, PA, USA) was subsequently infused (6 mg/kg) intravenously for 5 min, followed by 30 min stabilization time [4]. The blue (~440 nm) excitation of Oxyphor G2 allows assessment of microvascular pO2 in the renal cortex, in which it has been shown that hypoxic areas might arise during endotoxemia [11, 12].

The applied phosphorimeter and phosphorescent dye have been described elsewhere [10, 14]. Briefly, phosphorimetric sampling rate was set to 1 MHz (1 sample/μs). k 0 and k q were 1/251 and 273 μs−1 mmHg−1, respectively, corresponding to the constants previously found for Oxyphor G2 [10]. Phosphorimetric measurement time was set to 2/k 0, i.e., 502 μs.

Endotoxemic shock was induced by infusion of a bolus of lipopolysaccharide (LPS, 10 mg/kg, serotype 0127:B8, Sigma-Aldrich, St. Louis, MO, USA). When mean arterial pressure decreased to 40 mmHg, fluid resuscitation was started using a hydroxyethyl starch (HES) solution (5 ml/kg followed by 5 ml/kg/h; Voluven, 6% HES 130/0.4; Fresenius Kabi, Schelle, Belgium). During the entire protocol, heart rate, mean arterial pressure, end-tidal pCO2, and rectal temperature were monitored continuously. Microvascular pO2 was assessed every 10 s in the same renal phosphorimetric measurement volume (thus no spatial information) and microvascular pO2 histograms were recovered at baseline, during endotoxemic shock (LPS), and after 2 h of fluid resuscitation (LPS + HES).

3 Results

In the first set of simulations, we tested whether the application of a median filter or heterogeneity factor to the ESM would improve the recovery of pO2 histograms from noisy decay traces, having an SNR of 100 and 10. Simulated decay traces and corresponding histograms at a mean ± SD pO2 of 30 ± 15 and 90 ± 15 mmHg at SNRs of 100 and 10 are shown in Figs. 2, 3, 4. Noise does not affect the recovery of pO2 histograms with a low mean pO2 and widens the pO2 histograms with a high mean pO2, which could not be prevented with a median filter. All three methods adequately recover the input pO2 distributions both at high and low SNRs, i.e., input and recovered histograms are not significantly different (P > 0.05) as shown in Table 1. Overall, application of a median filter (Fig. 3) or heterogeneity factor (Fig. 4) to the ESM (Fig. 2) does not seem to improve the recovery of the histograms. Therefore, for future histogram analysis, no filter or heterogeneity factor will be applied.

Fig. 2
figure 2

Decay traces and corresponding histograms at a mean ± SD oxygen tension (pO2) of 30 ± 15 and 90 ± 15 mmHg at SNRs of 100 and 10. The input pO2 distributions are depicted by the black dashed line. The histograms are recovered using the exponential series method

Fig. 3
figure 3

Decay traces and corresponding histograms at a mean ± SD oxygen tension (pO2) of 30 ± 15 and 90 ± 15 mmHg at SNRs of 100 and 10. The input pO2 distributions are depicted by the black dashed line. The histograms are recovered using the exponential series method (ESM) with filtered decay traces

Fig. 4
figure 4

Decay traces and corresponding histograms at a mean ± SD oxygen tension (pO2) of 30 ± 15 and 90 ± 15 mmHg at SNRs of 100 and 10. The input pO2 distributions are depicted by the black dashed line. The histograms are recovered using the exponential series method (ESM) with the heterogeneity factor

Table 1 Kolmogorov-Smirnov test results (D values) comparing input and recovered oxygen tension (pO2) histograms

In the second set of simulations, we applied the ESM for the recovery of bimodal pO2 distributions, with a hypoxic and a normoxic mode at different ratios, at several SNRs and tested the ability of the analysis procedure to detect the extent of hypoxia. In Fig. 5, it is shown that, at high SNRs, the multi-exponential curve fitting analysis procedure adequately recovers bimodal distributions. As the noise level increases, the distinction between the two modes diminishes. However, the fraction of the first histogram bin remains linearly related to the degree of hypoxia (Fig. 5d), even at low SNRs.

Fig. 5
figure 5

Bimodal histograms with a hypoxic (0 ± 15 mmHg) and a normoxic (60 ± 15 mmHg) mode (ac). The contribution of the hypoxic mode was 10, 30, or 50% with signal-to-noise ratios (SNRs) of 50 (a), 20 (b), or 10 (c). In d, it is shown that the fraction of the first histogram bin is linearly related to the degree of hypoxia

To demonstrate the (patho)physiological relevance of multi-exponential curve fitting analysis of in vivo obtained phosphorescence decay traces for recovering microvascular pO2 histograms, measurements were performed in the rat kidney before and during endotoxemia. A typical baseline decay trace, having a SNR of 18 (=1/0.056), is depicted in Fig. 6a and the microvascular pO2 histograms at baseline, after induction of endotoxemia (LPS), and 2 h after the start of fluid resuscitation (LPS + HES) are shown in Fig. 6b. Mean microvascular oxygen tension, 〈pO2〉, obtained from the histograms was 93 ± 14 mmHg at baseline and decreased to 67 ± 13 mmHg after induction of endotoxemia. After 2 h of fluid resuscitation, 〈pO2〉 remained low at 64 ± 13 mmHg, but the hypoxic peak that had arisen during endotoxemia was reduced. Using mono-exponential curve fitting analysis the following values were obtained for microvascular pO2: 82 ± 19 mmHg at baseline, 46 ± 17 mmHg after induction of endotoxemia, and 47 ± 8 after 2 h of fluid resuscitation. Figure 6b, in combination with the results obtained by mono-exponential curve fitting analysis, illustrates the importance of recovering pO2 histograms under (patho)physiological conditions. The redistribution of oxygen by the treatment of endotoxemia with fluid resuscitation can be monitored and potentially provides insight into the nature of the renal insult caused by endotoxemia and the efficacy of its treatment with fluid resuscitation.

Fig. 6
figure 6

a A typical baseline decay trace, having a signal-to-noise ratio of 18 (=1/0.056). b Microvascular oxygen tension histograms at baseline, after induction of endotoxemia (LPS), and 2 h after the start of fluid resuscitation (LPS + HES). Mean microvascular oxygen tension, 〈pO2〉, was 93 ± 14 mmHg at baseline and decreased to 67 ± 13 mmHg after induction of endotoxemia. After 2 h of fluid resuscitation, 〈pO2〉 remained low at 64 ± 13 mmHg, but the hypoxic peak that had arisen during endotoxemia was reduced

4 Discussion and conclusions

The aim of this study was to extensively evaluate the ability of multi-exponential fitting analysis to adequately determine pO2 histograms from phosphorescence decay traces at different SNRs. For this purpose, we simulated decay traces representing uni- and bimodal pO2 distributions and recovered the pO2 histograms using (1) the ESM, (2) the ESM with filtered decay trace, and (3) the ESM with heterogeneity factor at several SNRs. Application of a median filter or heterogeneity factor to the ESM did not improve the recovery of the pO2 histograms. Ultimately, we applied the ESM for the recovery of microvascular pO2 histograms measured in the rat kidney in a model of endotoxemic shock and fluid resuscitation to demonstrate the (patho)physiological relevance of multi-exponential curve fitting analysis of in vivo obtained phosphorescence decay traces. We showed that the mean microvascular oxygen tension, 〈pO2〉, decreased after induction of endotoxemia and that after 2 h of fluid resuscitation, 〈pO2〉 remained low, but the hypoxic peak that had arisen during endotoxemia was reduced. This finding illustrates the importance of recovering pO2 histograms under (patho)physiological conditions as the redistribution of oxygen by the treatment of endotoxemia with fluid resuscitation can be monitored and potentially provides insight into the nature of the renal insult caused by endotoxemia and the efficacy of its treatment with fluid resuscitation [1113].

The (patho)physiological relevance of tissue pO2 heterogeneities has been appreciated for decades [36]. Oxygen electrodes have been extensively used to measure tissue pO2 at different sites and depths within organs to recover pO2 profiles [1, 17, 18, 22]. However, since oxygen electrode measurements are confined to a small tissue volume, only very superficial tissue oxygenation can be assessed noninvasively [32]. Mapping tissue oxygenation in depth requires insertion of one or more electrodes into the organ, causing cell damage and bleeding, which severely limits the applicability of this technique. The introduction of phosphorimetry [29, 30], in contrast, has opened the field of non-invasive and online monitoring of microcirculatory oxygenation and has initiated numerous studies investigating the (patho)physiological mechanisms underlying several types of shock [8] and the effects of different treatment strategies [e.g., 1113].

It is well established that endotoxemia-induced microvascular oxygenation heterogeneities have a key role in the development of numerous clinical complications, including renal failure [7, 1114]. However, although it is generally accepted that the oxygen-quenched phosphorescence and delayed fluorescence decay traces can be analyzed using multi-exponential curve fitting, only a few (patho)physiological studies have applied this methodology for the assessment of these critical oxygenation heterogeneities. This is probably because the reliability of the recovered pO2 histograms had never been extensively evaluated and lacked documentation. As the golden standard for measuring pO2 in vivo (i.e., the oxygen electrode) is unable to measure the pO2 heterogeneity, proper validation of multi-exponential curve fitting analysis, by comparison to this golden standard, is impossible. Hence, evaluating the reliability of this method can only be done using simulations.

With the introduction of the heterogeneity factor for the recovery of pO2 histograms by Golub et al. [6], the method was applied for recovering one Gaussian pO2 distribution at several SNRs. However, the ability of the ESM (with and without the heterogeneity factor) to recover pO2 histograms under noisy conditions depends on the pO2 distribution; low-pO2 distributions can be recovered adequately even under low SNRs, but high-pO2 distributions require relatively high SNRs. As the phosphorescence decay rate constant k (μs−1) is linearly related to the pO2, the number of data points above noise level, and thus fitting residuals, is inversely related to the pO2 and since the Levenberg–Marquardt procedure minimizes the sum of the squared residuals χ2 over the total extent of the decay trace. Thus, the relative contribution of the higher pO2 values to the fit is less than that of the lower pO2 values and therefore, noise affects the recovery of high pO2 values more than that of low pO2 values. A more extensive evaluation of this method, with both low and high pO2 distributions, was therefore required for determining its reliability.

Mik et al. [19] have performed such an evaluation of the ESM with heterogeneity factor and found that this method is indeed capable of recovering uni- and bimodal pO2 distributions for several SNRs. However, in that study the authors applied the heterogeneity factor as described by Golub et al. [6] and thereby amplified the noise at the end of the decay traces. Here, we present a modified version of the heterogeneity correction as developed by Golub and Zheng et al. and applied by Mik and Johannes et al. [6, 19, 20, 34]. Rather than transforming the phosphorescence decay trace using the heterogeneity factor, amplifying the noise in the latter part of the trace, we included the heterogeneity factor in the fit model. However, even though the noise was not amplified, in practice, the heterogeneity correction did not improve the recovery of pO2 histograms and was therefore omitted in the subsequent analyses. Furthermore, in attempt to reduce the influence of noise on the recovery of pO2 histograms, a median filter was applied. It should be noted, though, that application of a mean or median filter on noisy non-linear decay traces slightly ‘stretches’ the traces and thereby lowers the corresponding decay rates and thus pO2 values. As this ‘stretching’ property is stronger for the mean filter than for the median filter, we applied the latter (although this type of filter is usually applied to reduce the influence of outliers). However, application the median filter did not improve nor affect the recovery of pO2 histograms.

Some alternative methods are available for the recovery of pO2 histograms from phosphorescence and delayed fluorescence decay traces, each having specific advantages and disadvantages. The MEM, for instance, allows the recovery of very detailed histograms (typically 200 bins), but requires very high SNRs for stable operation (typically > 100) [16]. In our phosphorimetry setup, however, the measured phosphorescence decay traces typically have a SNR of ~20. Both the study of Golub et al. and that of Mik et al., have shown that the ESM operates stable at modest SNRs (i.e., ~20), but at the expense of reduced histogram detail [6, 19]. Recently, the MEM was questioned by Tsai et al. since the recovered pO2 distributions contained nonphysiologically high pO2 values, but this was later attributed to the reduced accuracy of the MEM for high pO2 values due to noise [27, 33]. Furthermore, in a study by Siemiarczuk et al., the ESM was compared to the MEM and was found to be significantly faster than the MEM [23]. The high temporal resolution of the ESM (in the order of seconds) is clearly advantageous as it allows rapid assessment of pO2 histograms online.

Recently, multi-exponential curve fitting analysis of oxygen-quenched phosphorescence decay traces, in combination with the Stern–Volmer equation, has been applied in several studies investigating microvascular and interstitial pO2 distributions in healthy and diseased tissue. Johannes et al. were able to identify the presence of ‘occult’ hypoxic microvascular areas during shock that would not have been detected when applying mono-exponential curve fitting analysis [11, 12]. This study supports these findings by showing that the detection of hypoxia with the ESM is highly adequate, even for very low SNRs. Golub and Pittman, moreover, studied cross-sectional pO2 profiles in microvasculature and found a bimodal distribution that was interpreted to arise from two separate physiological compartments within the phosphorimetric measurement volume, i.e., the intra- and extravascular space [5]. Using similar methodology, Mik et al. have explored pO2 heterogeneities in liver and heart mitochondria using multi-exponential curve fitting analysis of oxygen-quenched delayed fluorescence of intramitochondrial protoporphyrin IX [19, 20]. These studies underscore the (patho)physiological relevance of pO2 histograms and support their use in future research.

In conclusion, this study has characterized how noise affects the recovery of pO2 histograms using the ESM and documented the reliability of the ESM for recovering both low- and high-pO2 distributions for SNRs typically found in the experimental setting. This study might therefore serve as a frame of reference for investigations focused on oxygen (re)distribution during health and disease and encourage researchers to (re-)analyze data obtained in (earlier) studies possibly revealing new insights into complex disease states and treatment strategies. Moreover, now that the reliability of the ESM has been documented, new investigations focused on studying pO2 distributions under numerous (patho)physiological conditions could be initiated.