Skip to main content

Development of a mathematical model to estimate intra-tumor oxygen concentrations through multi-parametric imaging



Tumor hypoxia is involved in every stage of solid tumor development: formation, progression, metastasis, and apoptosis. Two types of hypoxia exist in tumors—chronic hypoxia and acute hypoxia. Recent studies indicate that the regional hypoxia kinetics is closely linked to metastasis and therapeutic responses, but regional hypoxia kinetics is hard to measure. We propose a novel approach to determine the local pO2 by fusing the parameters obtained from in vivo functional imaging through the use of a modified multivariate Krogh model.


To test our idea and its potential to translate into an in vivo setting through the use of existing imaging techniques, simulation studies were performed comparing the local partial oxygen pressure (pO2) from the proposed multivariate image fusion model to the referenced pO2 derived by Green’s function, which considers the contribution from every vessel segment of an entire three-dimensional tumor vasculature to profile tumor oxygen with high spatial resolution.


pO2 derived from our fusion approach were close to the referenced pO2 with regression slope near 1.0 and an r2 higher than 0.8 if the voxel size (or the spatial resolution set by functional imaging modality) was less than 200 μm. The simulation also showed that the metabolic rate, blood perfusion, and hemoglobin concentration were dominant factors in tissue oxygenation. The impact of the measurement error of functional imaging to the pO2 precision and accuracy was simulated. A Gaussian error function with FWHM equal to 20 % of blood perfusion or fractional vascular volume measurement contributed to average 7 % statistical error in pO2.


The simulation results indicate that the fusion of multiple parametric maps through the biophysically derived mathematical models can monitor the intra-tumor spatial variations of hypoxia in tumors with existing imaging methods, and the potential to further investigate different forms of hypoxia, such as chronic and acute hypoxia, in response to cancer therapies.


The oxygen concentration in normal tissue is maintained in a steady state, where the oxygen diffused from the capillaries meets the metabolic demand of the surrounding cells. Disrupting this balance, such as in wound repair for extra oxygen demand, triggers physiological angiogenesis. The nascent vessel sprouting and pruning from the existing vessel network in physiological angiogenesis are tightly regulated. If the disruption is severe enough, necrosis or apoptosis develops, resulting in a stroke or myocardial infarction. Unlike physiological angiogenesis, tumor-induced angiogenesis is uncontrolled, resulting in a tortuous, leaky, and dilated microvasculature. This dysfunctional microvasculature cannot supply the oxygen demand from tumor cells, and leads to a spatiotemporal heterogeneous hypoxia distribution in tumor [1]. The hypoxia further fuels abnormal angiogenesis and disrupts other cellular mechanisms, such as immune response [2] and cell survival [3, 4]. Thus, tumor hypoxia imaging and quantification are important in cancer research and have been an active research to search better cancer therapeutic outcomes [58]. There are two types of tumor hypoxia: chronic hypoxia (also known as diffusion-limited hypoxia) and acute hypoxia (also known as perfusion-limited hypoxia). Recent studies demonstrate that local tumor hypoxia behavior has profound impact in tumor metastasis and therapeutic outcomes [9, 10], hence characterizing and profiling the hypoxia in high spatiotemporal resolution can not only explore the further understanding causality between hypoxia and tumors but also is beneficial to cancer patient care.

Current in vivo preclinical diagnostic methods to determine the partial pressure of oxygen (pO2) include polarographic electrode (or Eppendorf probe), photoluminescence-quenching optical probe or biomarker assays based on the 2-nitroimidazoles (pimonidazole and EF5). These techniques provide a direct and quantifiable measure of pO2 or the relative oxygen level (pO2 < 10 mmHg) or hypoxic fraction in tissue, but their disadvantages have limited their wide spread use. For example, the Eppendorf probe necessitates an invasive procedure (needle insertion or biopsy) and lacks spatial information. Oxygen level quantification using 2-mitroimidazoles result in inconsistencies with Eppendorf probe measurements [11, 12] or with locoregional tumor control, disease-free survival [13] and event-free survival time [14]. PET and SPECT radiopharmaceutical tracers engineered around these biomarkers, such as 18F-MISO or 64Cu-ATSM, provide a noninvasive means to measure hypoxia [15, 16], and have been found to correlate with therapeutic response and outcome from radiation and chemotherapy [17, 18] when comparing average tumor-to-muscle ratios. However, these methods are still in the early stages of research and lack any mechanistic information causal to hypoxia.

Other approaches first quantify the hemodynamic parameters then extrapolate the local oxygen concentration. Near-infrared spectroscopy, such as diffuse optical tomography and photoacoustics, is widely used to measure hemoglobin concentration (CtHb) and its oxygen saturation (SaO2) level by quantifying the differential absorption spectrum of oxy- and deoxy-hemoglobin molecules [19, 20]. Quantitative BOLD MRI is also capable of measuring oxygen saturation levels in vivo as demonstrated in preclinical studies [21] and have been shown to correlate to pimonidizole staining. The oxygen concentration using dynamic contrast enhanced imaging by injecting contrast agents in blood circulation to measure local vascular hemodynamics correlates to Eppendorf pO2 histograph and pimonidizole-based immunohistochemistry [22, 23]. Each of these techniques provides a quantifiable measure of a single factor contributing to hypoxia, such as perfusion, hemoglobin concentration, or oxygen saturation, with good spatial and temporal resolution. By combining or fusing these parameters through the implementation of oxygen transport models, the precise local pO2 and etiology leading to hypoxia can be investigated. Integration of mathematical modeling and imaging observations proves to be a powerful approach to understand tumor vasculature behaviors at various scales and can have application in clinical diagnosis [2426]. The approach will become more reliable as the progression of imaging hardware provides accurate anatomical and functional measurements in high spatiotemporal resolution.

The advent of mathematical modeling in oxygen transport began as early as the twentieth century. Krogh in collaboration with Erlang, a mathematician, first approached the modeling of oxygen transport as a frame work to design experiments, interpret data, and suggest new insights on how the smallest micro vessels can supply oxygen to tissue, given the lack of technology at the time in support of his quest to understand oxygen transport in capillaries. With his model, he was able to conclude that the skeletal muscle tissue milieu was already well oxygenated and that the capillary bed only exchanged very low amount of oxygen with the muscle cells, contrary to the common belief at the time [27]. Krogh’s tissue cylindrical model set the stage for almost all subsequent models and has had tremendous and continuous impact to the field. By the 1960s, an interest in tissue oxygen transport was renewed, and steady-state and time-dependent models were developed to investigate the influence of hemoglobin and myoglobin oxygen binding within the microcirculation and whole in organs systems [28]. These models expanded upon Krogh’s model to include compartments for the RBCs, the plasma layer, the capillary wall, intestinal space and cellular volume, and used to study muscle physiology [29]. For example, Hellums used these concepts of mass transport to predict the passage of RBCs through capillaries would result in pO2 fluctuations [30], a phenomenon only recently observed [31, 32]. Multiple parallel vessel models better accounted for the O2 diffusion between capillaries, where advanced numerical techniques were beginning to applied, including finite difference and finite element methods [33]. When combined with advanced computational systems, local oxygen concentrations due to tissue heterogeneity could be further investigated [33]. By 2000, Secomb, Hsu, and Pries developed a Green’s function algorithm to calculate the steady-state local pO2 contributed from all vessel segments within a 3D network [34, 35], and was used to study the development of microvascular networks as well as different regulatory mechanisms on blood flow in tumors [3638]. By implementing finite difference methods, time-dependent solutions of heterogeneous microvascular networks under various physiological conditions were now possible [39].

The objective of this study is to profile the 3-D oxygen concentration and hypoxic fraction by fusing the multiple parameters obtained from in vivo functional imaging through the use of oxygen transport model. To obtain local oxygen concentration levels based on 3-D imaging, the Krogh model has been reformulated, where the inputs are obtained from the parametric images of the vascular physiology and hemoglobin status within the tumor. Thus, a multivariate image-fusion model of pO2 (or MVIF model) is provided. The MVIF model depends on the spatial resolution of the imaging system, as the point spread function decreases, the voxel volume will contain a single vessel and the local pO2 level will be determined with higher accuracy and precision, see Fig. 1. Prior to in vivo testing, simulation studies are performed to investigate the viability of the MVIF model, where previously published biological models and microscopic data were used as a point of reference. The goals of this study are to identify an optimal voxel size and the necessary instrumental sensitivities to obtain an accurate and precise determination of pO2. The range of spatial resolutions and sensitivities used in these simulation studies are consistent with in vivo dynamic contrast-enhanced (DCE) imaging (clinical and preclinical MRI, CT, or photoacoustic) and optical or photoacoustic spectroscopy (PCT-S). 3D tumor vasculature and hemodynamic quantification using intravital microscopy are used to simulate the in vivo parametric measurements (such as perfusion, fractional plasma volume, and the hemoglobin status within each voxel) which are subsequently fused to our model to obtain 3D pO2 maps. These MPO2 maps are compared on a voxel-by-voxel basis to the referenced pO2 values calculated by Green’s function [34]—which shows close correlation to the local tumor oxygen level using microelectrode [40]—to investigate the necessary image resolution and to determine the sensitivity by varying the statistical and systematic uncertainty in the in vivo imaging parameters.

Fig. 1
figure 1

Schematic illustration of the vascular reduction as function of the voxel size. As the voxel size is reduced, the vascularity of the regional vessel network within the voxel becomes less complex; thus, using a single effective vessel with multi vascular attributes to calculate the local pO2 becomes possible. In the MPO2 model, a single straight cylindrical vessel is placed at the center of the voxel and the average voxel pO2 (MPO2) is calculated. The effective vascular structure and its functional inputs (e.g. blood flow, hemoglobin concentration, and fraction of vessel volume) are obtained from in vivo imaging measurements, thus providing a means to fuse this information based on the biophysics model approach


Simulation procedure overview

The performance of our model was tested in two different tissue types: the brain and a tumor. The microvasculature within the brain served as an internal control or normal tissue, while the microvasculature of the tumor represented an abnormal or diseased tissue. The 3-D microvasculature of the tumor and the brain is displayed in Fig. 2a and b, the 3D dimension of tumor encompassing entire vasculature was 0.990 × 0.810 × 0.15 mm3, the dimension of brain was 0.150 × 0.160 × 0.14 mm3. Basically a window chamber was employed into the dorsal flap of the animal and tumor cells were implanted subcutaneously near a vessel in the chamber to study tumor-induced angiogenesis. Intravital microscopy was used to image 3D tumor vascular structure; photometric techniques were used to quantify hemodynamic variables (such as blood flow, hematocrit) in each vessel segment of the vasculature. The detailed experimental protocol and the procedures of post image and video processing to derive blood flow and vascular branching angle can be found in the publication of Fontanella et al. [41] and Brizel et al. [42]. The Green’s function algorithms to utilize the experimental vascular inputs from window chamber model in calculating oxygen distribution were downloaded from the website [43].

Fig. 2
figure 2

Illustration of reconstructed three-dimensional microvessel network of a a tumor (1000 × 1000× 150 μm3) and b brain (150 × 160 × 140 μm3) by confocal scanning microscope with sub-μm resolution (all axis’s are in units of micrometers)

In this study, the oxygen profile calculated by Green’s function was used as a reference. The Green’s function method calculates the pO2 at any location within the tumor by considering the contribution from all single vessel segments in the tissue volume. Once the pO2 profile within the entire tissue volume (brain or tumor) was completed, a virtual 3-D voxel grid representing the reconstructed image volume by DCE-CT was superimposed onto the tissue. The average effective hemodynamic values in each finite voxel were calculated based on its known microvascular architecture. These voxel-based hemodynamic parameters are used as inputs to our single-vessel multivariate oxygen transportation model (MVIF) to calculate the pO2 (MPO2) in each voxel of the virtual grid. The corresponding pO2 (GPO2) calculated using Green’s function was also calculated and used as a reference value. A scatter plot of MPO2 versus GPO2 for each voxel was compared by performing a linear regression. Since the placement of the virtual (or imaging) voxel grid is somewhat arbitrary and depends on the voxel size, the grid was randomly translated by less than a half a voxel length for larger voxel sizes. To account for these variations and to increase the number of data points for large-voxel MPO2 to GPO2 data, the grid was translated, the new vascular architecture for each voxel extracted, and the MPO2 and GPO2 calculated and plotted.

Implementing Green’s function algorithm to calculate referened pO2 (GPO2)

The Green’s function algorithm developed by Hsu and Secomb [34] was used to calculate the tissue oxygen concentrations resulting from the unique microvasculature geometry and hemodynamics in each tissue. The physiological constants used in these simulations were listed in Table 1. To reduce computation time, the tissue space was first discretized into an isotropic cubic lattice, the distance between two lattice points was 15 μm in our simulation (see Fig. 3). The oxygen consumption rate at each lattice point was derived from the Michaelis–Menten equation. Similarly vessels as oxygen source were discretized into sub vessel segments with 50 μm in length. The midpoint of each sub vessel segment represented the location oxygen source. In Green’s function algorithm, the major oxygen release factors (blood flow, oxygen diffusion rate in vessel and in tissue space, hemoglobin concentration, and oxygen release rate from hemoglobin into plasma simulated by Hill’s equation) and oxygen consumption rate were considered to calculate pO2 at the tissue lattice point (LPO2). Although Green’s function algorithm assumes the system to be under steady-state conditions and other assumptions to reduce the computational complexity, a statistically significant correlation between the tumor oxygen consumption rate derived by Green’s function and the actual microelectrode measurements in microvasculatures was reported [40], which provides a measure of confidence in pO2 from Green’s function as a reference.

Table 1 Key simulation parameters used in Green function model
Fig. 3
figure 3

A 2-D schematic of a cubic lattice used to evaluate the detailed pO2 profile in a tumor and the virtual grid to simulate the voxels from in vivo functional imaging modality. a The cubic lattice structure (blue diamonds) is used to calculate the pO2 value using Green’s function algorithm, where the edge length of the lattice is 15 µm. A 200-µm virtual grid is superimposed over the tumor. The color of the vessel segment indicates the geometrical relation of the vessel segment to the voxel. Green vessel segments are located entirely within the voxel; yellow segments represent vessels that cross two voxels; and blue segments are located outside the grid. b An example of a 100-µm virtual grid of the same tumor is displayed for comparison

Once all LPO2 was determined, a virtual grid composed of isotropic finite-sized voxels was superimposed onto the entire tissue, for example Fig. 3a and b demonstrated the virtual grid of 100- and 200-μm voxel overlaid on tumor. The voxel pO2, referred to as GPO2, was calculated by averaging all the LPO2 located inside the voxel, and represented the reference to evaluate the voxel pO2 by our multivariate single-vessel model.

Implementing the multivariate image-fusion (MVIF) model to calculate voxel pO2(MPO2)

All MVIF algorithms were programmed in Matlab (Mathworks Inc.) A detailed derivation of the multivariate image fusion model of oxygen concentration (MVIF) is described in the Additional file 1. MVIF is a modified single cylindrical vessel model proposed by Krogh and Erlang [27, 44] used to calculate the pO2 in a voxel of tissue based as a function of several key physiological parameters related to oxygen supply and consumption. The final analytic equation of MVIF to calculate pO2 in cylindrical coordinates is:

$$pO_{2} (r,z) = pO_{2} (z = 0) - \frac{M \cdot H}{(1 + m) \cdot F}G_{2} (z;r_{c} ,r_{T} ) - \frac{M \cdot H}{2D}G_{1} (r;r_{c} ,r_{T} ).$$

F represents blood perfusion, \(r_{c}\) is the single vessel radius, \(r_{T}\) is the upper radial boundary in MVIF, \(pO_{2}\,(z = 0)\) is the initial pO2 at the entrance of the vessel, \(D\) is the oxygen diffusion constant, \(H\) is Henry’s constant, \(M\) is the oxygen consumption rate, \(m\) is the slope of the oxygen-hemoglobin dissociation curve (see equation 2 in the Additional file 1), \(G_{1}\) is the geometrical term or form factor depicting radial diffusion within the voxel, and \(G_{2}\) is the geometrical term or form factor of advection along the z of cylindrical vessel.

In this simulation study, since the vasculature geometry and hemodynamic measurement inside the voxel of the virtual grid are known, some inputs for MVIF—such as the effective vessel radius \(r_{c}\), blood perfusion \(F\), and the initial pO2—can be calculated precisely as they were to measure by functional imaging modalities. Hill’s equation was used twice to calculate the initial pO2 in each voxel in this simulation. The lumen pO2 in each vessel segment of the vasculature was calculated using Green’s function, the corresponding SaO2 in each vessel segment, therefore, can be derived using Hill’s equation. In simulation the vessel segment ID and volume contribution within each voxel were recorded, the voxel SaO2 from all vessel segments within the voxel was calculated using the following formula:\(\sum\limits_{i} {{\text{SaO}}_{2} (i)*{\text{vol}}(i)/V}\), \({\text{SaO}}_{2} (i)\) is the oxygen saturation of vessel segment i; \({\text{vol}}(i)\) is the volume of vessel segment i within the voxel; \(V\) is the total vascular volume of the voxel. We used Hill’s equation again to derive the initial pO2 as input for MVIF from this voxel SaO2. The parameters or constants appearing in MVIF and Hill’s equation were listed in Table 2. The pO2 at any point within the voxel vessel was calculated; the average pO2 from the pO2 within the voxel represented the voxel pO2 calculated by MVIF, referred to as MPO2. Scattering plot of MPO2 versus GPO2 (as reference) on a voxel-by-voxel basis was plotted to evaluate the correlation and to identify the voxel outlier of which MPO2 significantly deviated from GPO2.

Table 2 Key simulation parameters used in MVIF model

MVIF response to voxel size (image spatial resolution)

In general, the complexity of the vessel network depends on the voxel size and the tissue type, as illustrated in Figs. 1, 2 and 3. In functional imaging image spatial resolution is determined by the voxel size. To investigate the accuracy and precision of the MVIF model as function of the voxels sizes, linear regression analysis was used to evaluate the correlation between MPO2 and GPO2 at different voxel sizes ranging from 50 to 300 µm, in 50 µm increments. For the brain tissue, only 50 and 100 µm voxel sizes were investigated because of the size of the tissue volume, 150 µm × 160 µm × 140 µm.

Hypoxic fraction correlation

The hypoxic fraction (HF) is the fraction of the tumor volume that has a pO2 less than a defined threshold, which typically ranged between 2.5 and 10 mmHg. Tumor HF has been used in clinical diagnosis, and has been estimated by a variety of techniques including immunohistochemistry, computed tomography, or [18F] 2-fluoro-2-deoxy-glocose PET [23, 45]. Similarly, the relationship between the HF as determined from MPO2 and GPO2 models was investigated. First, the tumor HF for MPO2 and GPO2 was calculated for various thresholds, from 2.5 to 15 mmHg in 2.5 mm steps. The number of voxels exceeding the hypoxic threshold was counted and divided by the total number of voxels in tumor. Since the brain (or normal) tissue is well oxygenated, a pO2 threshold ranging from 16 to 22 mmHg with 2 mmHg intervals was used.

Sensitivity and error analysis in MVIF modeling

Since nonlinear MVIF takes the hemodynamic and structural inputs from in vivo functional imaging to calculate MPO2, it is important to know (1) its response to changes in the input parameters and (2) how the hemodynamic or structural measurement error from in vivo functional imaging experiment affects the accuracy and precision of MPO2.

To answer the first objective the referenced MPO2 was calculated with the input values listed in Table 2. We then calculated MPO2 when only one input was considered as variable while other inputs remained as constant. Four new inputs which are ±10 and ±20 % from the referenced input value listed in Table 2 were used to calculate the corresponding MPO2 responses. The percent change in MPO2 compared to the referenced MPO2 value was then calculated and plotted. In this study we tested the MVIF response to following inputs: perfusion, vessel radius which is proportional to the fractional plasma volume, concentration of hemoglobin (CtHb), oxygen saturation (SaO2), and maximum oxygen consumption rate (M0). The reference values for blood perfusion, CtHb, and the SaO2 were selected from DCE-CT and PCT-S measurements [4649], and the maximum oxygen metabolic rate was chosen from the literature [50].

The second objective aims to learn how much the potential measurement error of functional DCE-CT could impact on the MVIF pO2 prediction in accuracy and precision. For this purpose we chose the voxel size equal to 150 μm because the scattering plot of GPO2 vs. MPO2 showed good correlation (see Fig. 5). Since blood perfusion and the fraction of vessel volume are two key DCE-CT measurements to tissue oxygenation, we evaluated the measurement error influence of those two major hemodynamic parameters to MVIF model. A Gaussian error function was used with a full-width at half-maximum (FWHM) equal to 20 % of the blood perfusion to generate ten mock blood perfusions for each voxel, and then the corresponding ten MPO2 were calculated. The average of ten MPO2 and the standard error were calculated and plotted as a function of GPO2. Identical procedures were applied to the fractional vessel volume to investigate how the error in the fraction of vessel measurements from DCE-CT could impact the MVIF performance.

Nearest-neighbor algorithm to reduce the outliers at 50-μM voxel scattering plot

The 50 μm voxel size for the tumor resulted in a significant increase in the variance and deviation from a slope of 1.0 as determined from the regression analysis (see Fig. 5a). This was not observed in the normal brain microvasculature. Unlike the brain, as the voxel size decreased for the tumor, an increasing number of voxels did not include a blood vessel, included a very small section of a blood vessel, or had a neighboring voxel with a large blood vessel. To compensate for the influence of the 6-neighboring voxels on a voxel’s average pO2 value, the average MPO2 from these surrounding voxels was added to the MPO2 value from the voxel itself to improve the correlation at 50 μm, or equivalently, to incorporate (or approximate) the boundary conditions. Therefore, the MPO2 calculation in each voxel was replaced by the average of surrounding six voxel pO2.


Correlation between MPO2 And GPO2 as a function of voxel size

For the normal brain microvasculature, the scatter plots of the average Green’s function pO2 (GPO2) and multivariate image fusion model pO2 (MPO2) for voxel dimensions of 50 and 100 µm were plotted in Fig. 4. At 50-μm voxel, the (slope, r2) values using linear regression analysis were (1.0, 0.86), at 100-μm voxel, the (slope, r2) were (1.0, 0.97). For the abnormal tumor microvasculature, the scatter plot and linear regression analysis for the six different voxel sizes (50–300 μm in 50 μm increment) were plotted in Fig. 5a–f. The (slope, r2) using linear regression analysis at voxel sizes from 50 to 300 μm in 50 μm increment were (0.99, 0.71), (0.98, 0.80), (1.0, 0.83), (1.0, 0.75), (0.65, 0.62), and (0.3, 0.46). The corresponding Pearson correlation coefficients (R) were 0.84, 0.89, 0.86, 0.84, 0.72 and 0.67 respectively, and all correlation coefficients reached statistical significance (p < 0.05).

Fig. 4
figure 4

Correlation between MPO2 and GPO2 in the brain. a Displayed is the scatter plot and linear regression fit for a a 50 μm and b 100 μm voxel size. The regression results (slope, offset, and r2) are shown in the top left corner

Fig. 5
figure 5

Correlation between MPO2 and GPO2 in the tumor at various voxel sizes. At 50-µm voxel size (a), two outlier populations (highlighted by green and red dash ovals) were observed. The green oval indicates the voxels that were overestimated, MPO2 significantly greater than GPO2, and the red oval the voxels that were underestimated, MPO2 significantly lower than GPO2. The correlation and linear regression analysis for voxel sizes 100, 150, and 200 μm (see bd) are better compared to the other voxel sizes (see a, e, f)

In Fig. 5a–c, the r2 value for the 50 μm voxel size was lower than that of 100 and 150 μm, which can be explained by the influence of oxygen contribution from the vessels (arterioles or venules) in the neighboring voxels. The nearest neighbor algorithm was applied to approximate this effect, where the average pO2 from the surrounding voxels was added to the MPO2 value and accounts voxel lacking vasculature and outliers observed in the scatter plots. After applying this algorithm, the correlation between MPO2 and GPO2 at 50 μm significantly improved (compare of Figs. 5a, 6a). The r2 as function of voxel size were displayed in Fig. 6b.

Fig. 6
figure 6

Correlation between MPO2 and GPO2 in the tumor after applying the nearest-neighbor algorithm for the 50-µm voxel data. a The scatter plot and linear regression fit are displayed after applying the nearest neighbor algorithm (see Fig. 5a for comparison). b The distribution of r2-values as a function of voxel size. The blue stars are from the r2-values in Fig. 5 and the red star is after applying the nearest neighbor algorithm to the 50-μm voxel data

Hypoxic fraction correlation between MPO2 And GPO2 as function of voxel size

Hypoxia fraction (HF) of tumor is another common clinical relevance parameter to evaluate tumor oxygen level. In the brain, the pO2 threshold values were set to 16, 18, 20, and 22 mmHg; and in the tumor, the thresholds were set to 2.5–15 mmHg in 2.5 mmHg increment. Figure 7 demonstrated the normal brain HF correlation and regression analysis of GPO2 and MPO2 model; the (slope, r2) for voxel size of 50 and 100 μm were (0.92, 0.97) and (0.96, 0.98) respectively. Figure 8 showed the tumor HF correlation and linear regression of voxel size 50–300 μm in 50 μm increment; the respective (slope, r2) of linear regression for voxel size from 50 to 250 μm in 50 μm increment were (0.44, 0.97), (0.84, 0.98), (0.85, 0.99), (0.79, 0.94) and (0.93, 0.86).

Fig. 7
figure 7

Normal brain tissue hypoxic fraction (HF) correlation of MPO2 and GPO2 in voxel sizes at 50 and 100 µm plotted in subplot a and b. In each voxel size, four HFs (16–22 mmHg in 2 mmHg increment) is used for linear regression analysis. The linear regression is plotted in red dash line; the result of linear regression and the Rsq values are displayed in top left corner

Fig. 8
figure 8

Tumor hypoxic fraction (HF) correlation of MPO2 and GPO2 in voxel sizes from 50 to 300 µm in 50 µm increment is plotted in subplot af. In each voxel size, six HFs (2.5–15 mmHg) is used for linear regression analysis. The linear regression is plotted in red dash line; the result of linear regression and the Rsq values are displayed in top left corner

Sensitivity of MVIF To microvascular hemodynamic and structural inputs

The changes in the local pO2, either GPO2 or MPO2, due to systematic offsets in the vascular inputs are displayed in Fig. 9. Those vascular inputs were chosen based on the reported influence on tumor oxygenation [51, 52], and determined from DCE-CT or PCT-S measurements [4649]. These values were 10 and 20 % change in the blood vessel radius (thus, fractional vascular volume), blood perfusion, oxygen metabolic rate, oxygen saturation (SaO2), and hemoglobin concentration (CtHb) relative to the reference values (Table 1). The absolute fluctuation in the average MPO2 from the four vascular inputs was 9.7, 15.2, 16.0, 12.6, and 15.0 %, respectively. The oxygen consumption rate had the largest impact on tissue oxygenation, while the vessel radius (or the fractional vascular volume) had the lowest impact on tissue oxygenation.

Fig. 9
figure 9

Sensitivity of MPO2 to the vascular input arguments. To investigate the systematic effects of changing the vascular arguments on pO2 in the tumor voxels, the vessel radius (Rc), blood perfusion, metabolic oxygen consumption rate, oxygen saturation rate (SaO2), and hemoglobin concentration (CtHb), are varied by −20 %, −10, 10, and 20 % from their respective reference values (Table 1). The percent deviation in the MPO2 relative to the original value (listed in Table 1) is plotted

MVIF response to the measurement error of microvascular inputs

To investigate how the precision of the measurement using in vivo functional imaging modality (e.g. DCE-CT) influence the MPO2 accuracy and precision, a Gaussian error function was used to generate ten configurations of perturbed maps of blood perfusions fractional vessel volume and the MPO2 and GPO2 maps calculated. The full width at half maximum (FWHM) of the Gaussian error function was set to 20 % of the mean vascular input value. The variation in MPO2 due to the measurement error in perfusion for the 150 μm voxel size is shown in Fig. 10a, where the error bar represents the standard deviation from these ten simulations. Similarly, the MPO2 variation due to the measurement error in the fractional vessel volume is shown in Fig. 10b. The Coefficient of variation (CV) is used to quantify the dispersion of the data relative to the mean. An average CV of the MPO2 as a result of the measurement error in blood perfusion and the fraction of vessel volume was 6.4 and 7.0 %.

Fig. 10
figure 10

The error in MPO2 due to the measurement uncertainty in the vascular input arguments. Plotted is the pO2 using MPO2 model versus GPO2 on a voxel-by-voxel basis, using the 100 μm voxel data. a Displayed is the simulated fractional vessel volume including the added noise using a Gaussian function (20 % noise at FWHM), and (b) and for perfusion


Many preclinical and clinical cancer investigations suggest that two types of hypoxia occur heterogeneously within a solid tumor: acute hypoxia (cycling hypoxia) and chronic hypoxia (diffusion-limited hypoxia) [10, 53]. In several tumor model acute hypoxia promotes tumor angiogenesis and tumor metastasis [10, 54], hence determining the etiology of hypoxia can provide new insights into predicting therapeutic response and developing novel therapeutic protocols to enhance efficacy. Here we proposed a novel assay to delineate diffusion- and perfusion-based hypoxia based on insufficient vasculature related to (G1 and D, the geometric diffusion length and rate, respectively) and poorly blood perfused tissue related to (G2 and F, the geometry related to convective forces and perfusion) under a steady-state conditions (see equation  3 and 4 in Supplement Data). By considering the variation of these parameters as a function of time [e.g., Perfusion = Perfusion(t) and m = m(t)], acute and chronic hypoxia could be assessed, which could be tested relative to variations in angiogenic (VEGF, FGF, etc.) and inflammatory (IL6, histamine, prostaglandins, NO) activity or their corresponding therapies. As an initial step, an in vivo high-resolution oxygen profiling assay is being proposed to image the spatial variations in tumor hypoxia. This capability will further advance our understanding the role hypoxia has in solid tumors as well as in developing new cancer treatments, such as in anti-angiogenic therapy.

In this study, a method was proposed that fuses microvascular functional attributes from functional imaging modalities through a multivariate oxygen transport model to obtain local tumor oxygen levels was investigated. High-resolution DCE-CT (or DCE-PCT) and PCT-S can longitudinally characterize local vessel structure and hemodynamics as well as tumor metabolic status. Those measurements represent tumor oxygenation status thus in principle a biophysical oxygen transportation model, in this case a multivariate modified Krogh model for a finite volume, can take the real measurements as inputs to calculate local oxygen level, and furthermore to profile tumor oxygen status. The major challenge to verify the model prediction is to have simultaneous tumor functional measurements and corresponding oxygen level measurement in the finite volume. To the best of our knowledge, no such experiment or data has been conducted or published. To evaluate this idea and our multivariate model, we used the well characterized tumor or brain microvasculature and functions using intravital confocal microscopy, and Secomb’s 3D oxygen transportation model to calculate the pO2 distribution as a reference to evaluate the voxel pO2 from MVIF model. The tumor oxygen consumption rates derived from Secomb’s 3-D model agree closely with the oxygen electrode probe measurements [40]. A major limitation of this study is the limited photon penetration from confocal microscopy, thus the simulation results from this study would only represent the peripheral microvasculature or microvasculature grown in window chamber [55, 56], which could be different from the microvasculature in other region of tumor. A second constraint of this investigation was the lack of few key measurements, such as the oxygen consumption rate, thus in simulation they were either derived from appropriated biophysical equation or assigned as constant.

As depicted in Fig. 1, the voxel size is associated with the vascular complexity or average inter vessel spacing and can have a significant impact on the accuracy and precision of the multivariate image-fusion algorithm (MVIF) model. An ideal model would provide a scatter plot with a flat r2 and a slope of 1.0 for all voxel sizes. As shown in Fig. 5b, this is clearly not the case, where an optimal voxel size was observed for the tumor and potentially the brain. Visually, this appears to depend on the voxel size relative to the inter vessel spacing or vascular density (see Fig. 3). For voxel sizes ranging from 100 to 200 μm, the coefficient of determination (r2), reached a maximum and remained relatively flat. For relatively small voxel sizes (50 μm), the average oxygen diffusion length and inter vessel spacing can exceed the voxel’s dimensions. As a result, two outlier groups were introduced into the scatter plot of MPO2 versus GPO2 for the tumor: an overestimated group, of which MPO2 was exceedingly greater than GPO2, and the underestimated group. Since these outliers were not observed for the 50 µm voxel data in the normal brain tissue, it demonstrates the impact of the abnormal and heterogeneous structure and function of the tumor vasculature on these results. To correct for these outliers, an algorithm was developed to account for oxygen diffusion from neighboring voxels (Fig. 6). For relatively larger voxel sizes, the precision and accuracy of MPO2 became worse and the slope (thus sensitivity) decreased; thus, the single vessel model can no longer adequately represent the complexity or heterogeneity of the vasculature, and as a result, the uncertainty in the average pO2 in a voxel becomes unacceptable. Given that the inter vessel spacing ranged from ~12 to 125 μm (Fig. 3), consistent with other tumor models (K, L, M), and was acquired within the periphery or rim of the tumor (Fig. 2), the results from this study show that the proposed modified Krogh image fusion model is accurate and precise for voxel sizes 200 μm and below (r2 < 0.75 and slope < 0.98). Therefore, the spatial resolution needed will need to be 200 μm or less. Such resolution can be obtained by many small animal imaging systems (photoacoustic, micro-CT, micro-MRI).

The hypoxic fraction (HF) indicates the hypoxic status of the entire tumor, a parameter used in 2-nitroimidazole based immunohistochemistry. Figure 7 showed the HF correlation of MPO2 and the referenced GPO2 with thresholds at 16, 18, 20, and 22 mmHg, a close regression was observed which shows that MVIF can accurately and precisely predict HF as a result of the normal microvasculature. In Fig. 8, the tumor HF data with threshold from 2.5 to 15.0 mmHg with 2.5 mmHg increment were linearly distributed; the linear regression analysis showed that the r2 at all voxel sizes was significant except at 300 µm albeit the HF from MVIF was consistently lower than that from Green’s function, and the regression offset was higher than that of brain microvasculature. Abnormal tumor microvessel structure or hemodynamics or both contributed to the HF offset and accuracy.

The oxygen metabolic rate, the blood flow perfusion, and hemoglobin concentration were three dominant factors in tissue oxygen transportation, as displayed in Fig. 9. Although the direct hemoglobin measurements and oxygen metabolic rates were not available in this simulation investigation, we chose Hill’s equation and Michaelis–Menton equation to estimate the hemoglobin concentration and oxygen metabolic rate in MPO2 calculation. For future in vivo experiments monitoring tumor pO2 fluctuations as function of space and time, all three parameters should be measured to achieve accurate and precise MPO2 estimation using MVIF. PCT-S is capable to derive hemoglobin concentration and oxygen saturation in high spatiotemporal resolution, the oxygen consumption rate can be derived from the fraction of cell volume from DCE-CT or dynamic contrast enhanced photoacoustic tomography. The use of photoacoustic imaging alone can avoid coregistration with DCE-CT, and the concern of radiation superimposing on the tumor progression or compound treatment.

The influence of measurement error of functional imaging to the MPO2 precision and accuracy was simulated and summarized in Fig. 10. Overall, 20 % error (FWHM) in perfusion or fractional vessel volume did not deteriorate accuracy and precision when at least 10 measurements were acquired at each voxel and used to estimate the pO2.

There remains limitations of this technique as well as some potential future directions. A key limitation of the proposed technique is for larger voxel sizes. If a model can be devised to rescue the MPO2 versus GPO2 correlation beyond 200 μm, the ability to translate into the clinic is viable.

Multi-vessel models consisting of 2 and 3 parallel cylindrical vessels with identical dimensions and hemodynamics and uniformly spaced within 250 and 300 µm voxels were simulated. However, the improvement was limited. Introducing symmetry into these models and simple bifurcated structures could explain the scatter and nonlinearity in the data; however, identifying or imposing a single concept or methodology would challenging, and requires additional constraints, such as including angiogenic information or inferring a scaling factor based on fractals [57, 58]. Additional pathological features affecting oxygen transport should also be considered, such as oxygen permeability across the vessel wall. Even though tumor vasculature is highly fenestrated and the pressure gradient across the vessel wall reduced due to elevated interstitial fluid pressure, models including oxygen permeability suggest this may have a non-negligible effect [59, 60]. With the ability to measure the permeability-surface area product and interstitial fluid pressure through in vivo imaging, it may be possible to test these models [46, 61]. Green’s function solution of oxygen transport has a number of advantages in the way it handles boundary conditions, the ability to handle a larger set of physiological parameters, efficiency of calculation, and validation to in vivo measurements, it lacks the ability to handle necessary advanced features FDTD and FEM methods support, such as dynamic changes in its parameters and oxygen permeability.


The results from this simulation study demonstrate that the fusion of in vivo functional imaging based on the MPO2 model to quantify local tumor oxygen concentrations is feasible. A significant correlation was measured between the MPO2 model, a single vessel model, and the Green’s function algorithm (GPO2), a detailed microvascular model, for voxel sizes ranged from 50 to 200 μm, where MPO2 was calculated based on in vivo functional imaging measurements. This upper limit is consistent with the spatial resolution from existing small animal scanners, and provides a new technique to assay in vivo intra-tumor hypoxia. Future work will investigate the potential to monitor the dynamic changes in pO2. This additional capability will be able to characterize regions of the tumor undergoing chronic and acute forms of hypoxia and the hemodynamics parameters responsible, thus providing critical information on the role hypoxia and metabolism play in cancer progression and therapy, in particular anti-angiogenic therapies [10, 62].



dynamic contrast-enhanced computed tomography


photoacoustic computed tomographic spectroscopy

pO2 :

oxygen partial pressure in mmHg

CtHb :

hemoglobin concentration

SaO2 :

oxygen saturation


average pO2 derived in a finite volume using Green’s function method


multivariate image–fusion algorithm


average pO2 derived in a finite volume using MVIF model

Rsq or r2 :

R-squared value


  1. Jain RK. Normalization of tumor vasculature: an emerging concept in antiangiogenic therapy. Science. 2005;307(5706):58–62.

    Article  Google Scholar 

  2. Barsoum IB, et al. A mechanism of hypoxia-mediated escape from adaptive immunity in cancer cells. Cancer Res. 2014;74(3):665–74.

    Article  Google Scholar 

  3. Huang Y, et al. Vascular normalization as an emerging strategy to enhance cancer immunotherapy. Cancer Res. 2013;73(10):2943–8.

    Article  Google Scholar 

  4. Coussens LM, Zitvogel L, Palucka AK. Neutralizing tumor-promoting chronic inflammation: a magic bullet? Science. 2013;339(6117):286–91.

    Article  Google Scholar 

  5. Semenza GL. Hypoxia-inducible factors: mediators of cancer progression and targets for cancer therapy. Trends Pharmacol Sci. 2012;33(4):207–14.

    Article  Google Scholar 

  6. Hendrickson K, et al. Hypoxia imaging with [F-18] FMISO-PET in head and neck cancer: potential for guiding intensity modulated radiation therapy in overcoming hypoxia-induced treatment resistance. Radiother Oncol. 2011;101(3):369–75.

    Article  Google Scholar 

  7. Bell C, et al. Hypoxia imaging in gliomas with 18F-fluoromisonidazole PET: toward clinical translation. Semin Nucl Med. 2015;45(2):136–50.

    Article  Google Scholar 

  8. Cao N, et al. Monitoring the effects of anti-angiogenesis on the radiation sensitivity of pancreatic cancer xenografts using dynamic contrast-enhanced computed tomography. Int J Radiat Oncol Biol Phys. 2014;88(2):412–8.

    Article  Google Scholar 

  9. Wilson WR, Hay MP. Targeting hypoxia in cancer therapy. Nat Rev Cancer. 2011;11(6):393–410.

    Article  Google Scholar 

  10. Rofstad EK, et al. Tumors exposed to acute cyclic hypoxic stress show enhanced angiogenesis, perfusion and metastatic dissemination. Int J Cancer. 2010;127(7):1535–46.

    Article  Google Scholar 

  11. Nordsmark M, et al. Measurements of hypoxia using pimonidazole and polarographic oxygen-sensitive electrodes in human cervix carcinomas. Radiother Oncol. 2003;67(1):35–44.

    Article  Google Scholar 

  12. Evans SM, et al. Comparative measurements of hypoxia in human brain tumors using needle electrodes and EF5 binding. Cancer Res. 2004;64(5):1886–92.

    Article  Google Scholar 

  13. Kaanders JH, et al. Pimonidazole binding and tumor vascularity predict for treatment outcome in head and neck cancer. Cancer Res. 2002;62(23):7066–74.

    Google Scholar 

  14. Evans SM, et al. Patterns and levels of hypoxia in head and neck squamous cell carcinomas and their relationship to patient outcome. Int J Radiat Oncol Biol Phys. 2007;69(4):1024–31.

    Article  Google Scholar 

  15. Mees G, et al. Molecular imaging of hypoxia with radiolabelled agents. Eur J Nucl Med Mol Imaging. 2009;36(10):1674–86.

    Article  Google Scholar 

  16. Troost EG, et al. Correlation of [18F]FMISO autoradiography and pimonidazole [corrected] immunohistochemistry in human head and neck carcinoma xenografts. Eur J Nucl Med Mol Imaging. 2008;35(10):1803–11.

    Article  Google Scholar 

  17. Matsumoto K, et al. The influence of tumor oxygenation on hypoxia imaging in murine squamous cell carcinoma using [64Cu]Cu-ATSM or [18F]Fluoromisonidazole positron emission tomography. Int J Oncol. 2007;30(4):873–81.

    Google Scholar 

  18. Padhani AR, et al. Imaging oxygenation of human tumours. Eur Radiol. 2007;17(4):861–72.

    Article  Google Scholar 

  19. Yalowitz JA, et al. Cytotoxicity and cellular differentiation activity of methylenebis(phosphonate) analogs of tiazofurin and mycophenolic acid adenine dinucleotide in human cancer cell lines. Cancer Lett. 2002;181(1):31–8.

    Article  Google Scholar 

  20. Bayer CL, Luke GP, Emelianov SY. Photoacoustic imaging for medical diagnostics. Acoust Today. 2012;8(4):15–23.

    Article  Google Scholar 

  21. Christen T, et al. Tissue oxygen saturation mapping with magnetic resonance imaging. J Cereb Blood Flow Metab. 2014;34(9):1550–7.

    Article  Google Scholar 

  22. Cooper RA, et al. Tumour oxygenation levels correlate with dynamic contrast-enhanced magnetic resonance imaging parameters in carcinoma of the cervix. Radiother Oncol. 2000;57(1):53–9.

    Article  Google Scholar 

  23. Gulliksrud K, et al. Quantitative assessment of hypoxia in melanoma xenografts by dynamic contrast-enhanced magnetic resonance imaging: intradermal versus intramuscular tumors. Radiother Oncol. 2010;97(2):233–8.

    Article  Google Scholar 

  24. Kim E, et al. Multiscale imaging and computational modeling of blood flow in the tumor vasculature. Ann Biomed Eng. 2012;40(11):2425–41.

    Article  Google Scholar 

  25. Stamatelos SK, et al. A bioimage informatics based reconstruction of breast tumor microvasculature with computational blood flow predictions. Microvasc Res. 2014;91:8–21.

    Article  Google Scholar 

  26. Toma-Dasu I, Dasu A, Brahme A. Quantifying tumour hypoxia by PET imaging—a theoretical analysis. Adv Exp Med Biol. 2009;645:267–72.

    Article  Google Scholar 

  27. Krogh A. The number and distribution of capillaries in muscles with calculations of the oxygen pressure head necessary for supplying the tissue. J Physiol. 1919;52(6):409–15.

    Article  Google Scholar 

  28. Popel AS. Theory of oxygen transport to tissue. Crit Rev Biomed Eng. 1989;17(3):257–321.

    Google Scholar 

  29. Eggleton CD, et al. Calculations of intracapillary oxygen tension distributions in muscle. Math Biosci. 2000;167(2):123–43.

    Article  MATH  Google Scholar 

  30. Hellums JD. The resistance to oxygen transport in the capillaries relative to that in the surrounding tissue. Microvasc Res. 1977;13(1):131–6.

    Article  Google Scholar 

  31. Golub AS, Pittman RN. Erythrocyte-associated transients in PO2 revealed in capillaries of rat mesentery. Am J Physiol Heart Circ Physiol. 2005;288(6):H2735–43.

    Article  Google Scholar 

  32. Barker MC, Golub AS, Pittman RN. Erythrocyte-associated transients in capillary PO2: an isovolemic hemodilution study in the rat spinotrapezius muscle. Am J Physiol Heart Circ Physiol. 2007;292(5):H2540–9.

    Article  Google Scholar 

  33. Popel AS, Charny CK, Dvinsky AS. Effect of heterogeneous oxygen delivery on the oxygen distribution in skeletal muscle. Math Biosci. 1986;81(1):91–113.

    Article  MATH  Google Scholar 

  34. Secomb TW, et al. Green’s function methods for analysis of oxygen delivery to tissue by microvascular networks. Ann Biomed Eng. 2004;32(11):1519–29.

    Article  Google Scholar 

  35. Secomb TW, et al. Theoretical simulation of oxygen transport to brain by networks of microvessels: effects of oxygen supply and demand on tissue hypoxia. Microcirculation. 2000;7(4):237–47.

    Article  Google Scholar 

  36. Secomb TW, et al. Analysis of oxygen transport to tumor tissue by microvascular networks. Int J Radiat Oncol Biol Phys. 1993;25(3):481–9.

    Article  Google Scholar 

  37. Secomb TW. Theoretical models for regulation of blood flow. Microcirculation. 2008;15(8):765–75.

    Article  Google Scholar 

  38. Pries AR, et al. The shunt problem: control of functional shunting in normal and tumour vasculature. Nat Rev Cancer. 2010;10(8):587–93.

    Article  Google Scholar 

  39. Goldman D, Popel AS. A computational study of the effect of capillary network anastomoses and tortuosity on oxygen transport. J Theor Biol. 2000;206(2):181–94.

    Article  Google Scholar 

  40. Dewhirst MW, et al. Determination of local oxygen consumption rates in tumors. Cancer Res. 1994;54(13):3333–6.

    Google Scholar 

  41. Fontanella AN, et al. Quantitative mapping of hemodynamics in the lung, brain, and dorsal window chamber-grown tumors using a novel, automated algorithm. Microcirculation. 2013;20:724–35.

    Google Scholar 

  42. Brizel DM, et al. A comparison of tumor and normal tissue microvascular hematocrits and red cell fluxes in a rat window chamber model. Int J Radiat Oncol Biol Phys. 1993;25(2):269–76.

    Article  Google Scholar 

  43. Secomb TW. Accessed on May 26 2011.

  44. Krogh A. The rate of diffusion of gases through animal tissues, with some remarks on the coefficient of invasion. J Physiol. 1919;52(6):391–408.

    Article  Google Scholar 

  45. Komar G, et al. 18F-EF5: a new PET tracer for imaging hypoxia in head and neck cancer. J Nucl Med. 2008;49(12):1944–51.

    Article  Google Scholar 

  46. Cao M, et al. Developing DCE-CT to quantify intra-tumor heterogeneity in breast tumors with differing angiogenic phenotype. IEEE Trans Med Imaging. 2009;28(6):861–71.

    Article  Google Scholar 

  47. Krishnamurthi G, et al. Functional imaging in small animals using X-ray computed tomography–study of physiologic measurement reproducibility. IEEE Trans Med Imaging. 2005;24(7):832–43.

    Article  Google Scholar 

  48. Stantz KM, Liang Y, Hutchins GD. Kinematic modeling and its implication in longitudinal chemotherapy study of tumor physiology: ovarian xenograft mouse model and contrast-enhanced dynamic CT. Proc SPIE. 2004;5369:769–79.

    Article  Google Scholar 

  49. Stantz KM, et al. Photoacoustic spectroscopic imaging of intra-tumor heterogeneity and molecular identification. Proc SPIE BIOS. 2006;6086:36–47.

    Google Scholar 

  50. Arthur C, Guyton JEH. Textbook of medical physiology. Philadelphia: Elsevier Saunders; 2006.

    Google Scholar 

  51. Batchelor TT, et al. Improved tumor oxygenation and survival in glioblastoma patients who show increased blood perfusion after cediranib and chemoradiation. Proc Natl Acad Sci USA. 2013;110(47):19059–64.

    Article  Google Scholar 

  52. Stylianopoulos T, Jain RK. Combining two strategies to improve perfusion and drug delivery in solid tumors. Proc Natl Acad Sci USA. 2013;110(46):18632–7.

    Article  Google Scholar 

  53. Palmer GM, et al. Optical imaging of tumor hypoxia dynamics. J Biomed Opt. 2010;15(6):066021.

    Article  Google Scholar 

  54. Rofstad EK, et al. Fluctuating and diffusion-limited hypoxia in hypoxia-induced metastasis. Clin Cancer Res. 2007;13(7):1971–8.

    Article  Google Scholar 

  55. Palmer GM, et al. In vivo optical molecular imaging and analysis in mice using dorsal window chamber models applied to hypoxia, vasculature and fluorescent reporters. Nat Protoc. 2011;6(9):1355–66.

    Article  Google Scholar 

  56. Papenfuss HD, et al. A transparent access chamber for the rat dorsal skin fold. Microvasc Res. 1979;18(3):311–8.

    Article  Google Scholar 

  57. Gazit Y, et al. Fractal characteristics of tumor vascular architecture during tumor growth and regression. Microcirculation. 1997;4(4):395–402.

    Article  Google Scholar 

  58. Baish JW, Jain RK. Fractals and cancer. Cancer Res. 2000;60(14):3683–8.

    Google Scholar 

  59. Skeldon AC, et al. Modelling and detecting tumour oxygenation levels. PLoS ONE. 2012;7(6):e38597.

    Article  Google Scholar 

  60. Barrett MJ, Suresh V. Extra permeability is required to model dynamic oxygen measurements: evidence for functional recruitment? J Cereb Blood Flow Metab. 2013;33(9):1402–11.

    Article  Google Scholar 

  61. Hompland T, et al. Assessment of the interstitial fluid pressure of tumors by dynamic contrast-enhanced magnetic resonance imaging with contrast agents of different molecular weights. Acta Oncol. 2013;52(3):627–35.

    Article  Google Scholar 

  62. Cairns RA, Hill RP. Acute hypoxia enhances spontaneous lymph node metastasis in an orthotopic murine model of human cervical carcinoma. Cancer Res. 2004;64(6):2054–61.

    Article  Google Scholar 

  63. Grote J, Susskind R, Vaupel P. Oxygen diffusivity in tumor tissue (DS-carcinosarcoma) under temperature conditions within the range of 20–40 degrees C. Pflugers Arch. 1977;372(1):37–42.

    Article  Google Scholar 

  64. Zander R. Cellular oxygen concentration. Adv Exp Med Biol. 1975;75:463–7.

    Article  Google Scholar 

  65. Fournier RL. Basic transport phenomena in biomedical engineering. 3rd ed. Boca Raton: CRC Press; 2011.

    Google Scholar 

Download references

Authors’ contributions

KMS carried out the idea in fusing mathematical models and in vivo functional imaging to estimate entire 3D oxygen distribution in tumor; CWL carried out simulation programming and evaluation using oxygen profile from Green’s function as a standard. All authors read and approved the final manuscript.


The first author would like to give thanks to the support from previous employer: Eli Lilly and Company, and Covance. The first author also would like give thanks to the support from School of Health Sciences at Purdue University in West Lafayette IN.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

The datasets during and/or analysed during the current study available from the corresponding author on reasonable request.


NIH/NIBIB R21 EB012849, PI: Stantz KM.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Keith M. Stantz.

Additional file

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lee, CW., Stantz, K.M. Development of a mathematical model to estimate intra-tumor oxygen concentrations through multi-parametric imaging. BioMed Eng OnLine 15, 114 (2016).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: