## Abstract

Oxygen levels in cancerous tissue can have a significant effect on treatment response: hypoxic tissue is both more radioresistant and more chemoresistant than well-oxygenated tissue. While recent advances in medical imaging have facilitated real-time observation of macroscopic oxygenation, the underlying physics limits the resolution to the millimetre domain, whereas oxygen tension varies over a micrometre scale. If the distribution of oxygen in the tumour micro-environment can be accurately estimated, then the effect of potential dose escalation to these hypoxic regions could be better modelled, allowing more realistic simulation of biologically adaptive treatments. Reaction–diffusion models are commonly used for modelling oxygen dynamics, with a variety of functional forms assumed for the dependence of oxygen consumption rate (OCR) on cellular status and local oxygen availability. In this work, we examine reaction–diffusion models of oxygen consumption in spherically and cylindrically symmetric geometries. We consider two different descriptions of oxygen consumption: one in which the rate of consumption is constant and one in which it varies with oxygen tension in a hyperbolic manner. In each case, we derive analytic approximations to the steady-state oxygen distribution, which are shown to closely match the numerical solutions of the equations and accurately predict the extent to which oxygen can diffuse. The derived expressions relate the limit to which oxygen can diffuse into a tissue to the OCR of that tissue. We also demonstrate that differences between these functional forms are likely to be negligible within the range of literature estimates of the hyperbolic oxygen constant, suggesting that the constant consumption rate approximation suffices for modelling oxygen dynamics for most values of OCR. These approximations also allow the rapid identification of situations where hyperbolic consumption forms can result in significant differences from constant consumption rate models, and so can reduce the computational workload associated with numerical solutions, by estimating both the oxygen diffusion distances and resultant oxygen profile. Such analysis may be useful for parameter fitting in large imaging datasets and histological sections, and allows easy quantification of projected differences between functional forms of OCR.

## 2. Introduction

It has been known since the pioneering work of Gray and co-workers in the 1950s [1] that well-oxygenated tissue responds better to radiotherapy and chemotherapy than by up to a factor of 3 relative to tumours with extensive hypoxia. This relative boosting effect of oxygen on cell-kill is often quantified by the oxygen enhancement ratio (OER), which is of fundamental importance in radiotherapy [2]. The prognostic influence of hypoxia has led to the concept of dose painting [3,4] in radiotherapy, which proposes that hypoxic regions of a tumour could be given an increased dose to mitigate their inherent radioresistance. Imaging modalities such as positron emission tomography (PET) can be coupled with hypoxia binding tracers such as 18F-fluoromisonidazole to allow the non-invasive estimation of hypoxia *in vivo* [3]. Although this approach is promising, the spatial resolution remains prohibitively low [5], as PET imaging is physically constrained to a millimetre-scale resolution, while oxygen tension is known to vary significantly over tens of micrometres [6–9], owing to the relatively short oxygen diffusion distances in tissue.

To maximize the potential advantages of a dose painting approach to treatment outcome, it is vital that the underlying oxygen distributions are well understood. Mathematical modelling of oxygen transport at micrometre scales is therefore required to estimate oxygen distributions. In continuum formulations, oxygen is assumed to diffuse through tissue while simultaneously being consumed by the respiring cells. Models of this type have been developed for spherical geometries to describe the oxygen tension in multicellular tumour spheroid (MCTS) growth [10–12] and cylindrical geometries for blood vessel perfusion modelling [8,9,13,14]. This approach has been widely used in various tissue geometries and situations [8–11,15–20] to estimate the oxygen dynamics through a tissue. Several functional forms of the oxygen consumption rate (OCR) have been proposed. In some formulations, it is treated as a constant [8,9] while elsewhere the OCR is assumed to vary with oxygen tension, typically obeying a rectangular hyperbolic relationship, similar in form to Michaelis–Menten kinetics or the Langmuir adsorption isotherm. This particular functional form has been shown to be a good approximation in a variety of tissues, both healthy and cancerous [21–24]. A typical hyperbolic reaction rate is given by
*s* is the substrate concentration, *K*_{m} is substrate level where *ν* is half

For cellular oxygen consumption, this model is phenomenological, and its sigmoidal form is intended to capture the observation that the OCR in a given tissue is typically lower under hypoxic conditions than under normoxia. From a biological perspective, such dynamics may be explained by the observation that under hypoxic conditions, expression of proteins such as HIF-1 can lead to downregulation of mitochondrial oxygen consumption, as well as affecting cell cycling [25]. The nonlinear form of (2.1) renders reaction–diffusion equations involving hyperbolic kinetics analytically intractable in general, and so they must be solved numerically.

The form of equation (2.1) suggests that the relative magnitude of the constant *K*_{m} for oxygen has a distinct effect on the shape of the consumption rate curve. The literature to date indicates that this constant for oxygen is quite low compared with typical oxygen tension in a vessel, typically *K*_{m}≃1 mm Hg [26–28]. This may suggest only small predicted differences between the two consumption rate forms described above in most cases, though there is no simple method for estimating which combinations of parameters may result in a significant difference between both forms. More importantly perhaps, the oxygen diffusion limit could differ significantly between both models; consider two reaction–diffusion equations with a maximum OCR of *a*, the first operating under the assumption of constant OCR and the second under the assumption of hyperbolic kinetics. For the first model, OCR is intrinsic and is not modulated by oxygen availability so oxygen consumption is *a* throughout. In the second model, the OCR is modulated by the presence of oxygen and is given by *aP*/(*P*+*K*_{M}), where *P* is the local oxygen partial pressure. In this case, oxygen consumption will reduce with local oxygen availability, allowing oxygen to diffuse deeper into the respiring tissue than in the constant case. Thus, if all other parameters are kept the same, the oxygen diffusion distance in the hyperbolic model will be always greater than for the constant consumption rate case. Under what circumstances this difference may be significant in a given situation is an open question, and part of the motivation for this work.

Previous literature has outlined how oxygen diffusion distance can be analytically derived for the constant consumption rate case in the case of cylindrical [8] and spherical [10] geometries, and to date numerical approaches have been used when considering non-constant kinetics, owing to the analytic intractability of the problem [29]. The significant drawback of such approaches is that the diffusion boundaries must be specified by the user, rendering the process uncertain and making direct comparison between constant oxygen consumption models and hyperbolic kinetic models difficult. Despite the small literature estimates of *K*_{m} for oxygen, there may be cases where significant differences emerge between constant consumption and M–M case. It would be beneficial to have a simple method of quantifying the effects of both constant OCR and hyperbolic consumption rate forms, and a specifically a method for estimating the predicted oxygen diffusion distance in the hyperbolic case explicitly. Such a method would also have to be capable of explicitly relating OCR to the oxygen diffusion distance. In this work, we derive analytical approximations for hyperbolic-type models in both spherical and cylindrical geometry, and contrast these to constant consumption models. We show how oxygen diffusion distance is governed by OCR, and we show how the models derived can be used to obtain a value for the oxygen diffusion distance and rapid estimation of the expected oxygen profile. As the oxygen tension and associated boundaries can be estimated analytically, this can reduce computational workload when modelling systems of such elements or performing image analysis on vessel sections. We present results of this analysis for key quantities of interest such as expected diffusion distance and oxygen distribution and discuss the implications for oxygen modelling at a micrometre scale.

## 3. Methods and models

### 3.1 Spherical symmetry

We adopt a continuum modelling approach, deriving a reaction–diffusion equation governing the spatio-temporal dynamics of oxygen in an avascular tissue. Adopting spherical polar coordinates, we suppose that the centre of the MCTS is located at *r*=0. Since the average doubling time of cells in an MCTS is much longer than the typical time scale of oxygen diffusion across the spheroid, we make the simplifying assumption that the spheroid radius, denoted by *r*_{o}, is constant over the time scale of measurement.

We initially suppose that the spheroid is sufficiently large that a central anoxic region exists where the oxygen concentration is zero. We denote the radius of this anoxic region by *r*_{n}, so that oxygen partial pressure is zero in the region 0≤*r*≤*r*_{n}, and *r*_{n} is the boundary to which oxygen can penetrate from *r*_{o}. For convenience, we let *r*_{c}=*r*_{o}−*r*_{n} denote the width of the viable rim of the tissue. We further assume the MTCS is suspended in a well-mixed growth medium [10,30].

We assume that spherical symmetry is preserved, as illustrated in figure 1, and we let *p*(*r*,*t*) denote the partial pressure of oxygen at a distance *r* from the centre of the MCTS at time *t*. We consider changes in oxygen partial pressure due to two processes: Fickian diffusion into and within the MCTS, and cellular consumption. We also make the simplifying assumption that the diffusion coefficient of oxygen, denoted by *D*, is spatially invariant. We let *a*(*p*) denote the rate of cellular consumption of oxygen, which may depend on the local oxygen partial pressure. It can readily be shown that for typical literature values of the oxygen diffusion constant *D* that oxygen diffusion through an MCTS occurs on a time scale of seconds and a steady-state approximation is therefore acceptable, under which *p* satisfies the reaction–diffusion equation
*r*_{n}<*r*<*r*_{o}, where *Ω*=3.0318×10^{7} mm Hg kg m^{−3} arises from Henry’s law as outlined in previous work [10]. Although one could readily define a constant that subsumes *a*(*p*) and *Ω*, we will keep them separate in this work to reflect their physical meaning—specifically *a*(*p*) is the oxygen consumption rate in SI units of volume of oxygen per unit mass per unit time. If the two terms are subsumed into a resultant term *Ωa*(*p*), then this yields consumption in units of oxygen partial pressure per unit time. To close this system, we must impose appropriate boundary conditions. We impose the regularity condition *p*(*r*_{n})=0 at the anoxic edge and assume that the exterior medium is well mixed and that oxygen is in excess there. This is equivalent to assuming a very high mass transfer coefficient for *p* at *r*=*r*_{o} leading to the Dirichlet boundary condition *p*(*r*_{o})=*p*_{o}. It is important to note that *r*_{n}, the anoxic radius, is not known *a priori* but can be derived from first principles, as outlined below.

#### 3.1.1 Constant consumption rate

If the OCR is taken to be constant then it is possible to derive a closed-form analytic solution to the model [10]. In this case, oxygen consumption is a constant so that *a*(*p*)=*a*_{o}. The exact solution is given analytically by
*r*_{n}≤*r*≤*r*_{o}. It can further be shown [10] that the maximum radius that oxygen can penetrate before the formation of an anoxic core is given by
*r*_{n} is not known *a priori* and needs to be calculated. In equation (3.2), the anoxic radius *r*_{n} is undetermined and can be explicitly calculated [10] to give
*r*_{n}=0 and a similar analysis can be performed.

#### 3.1.2 Hyperbolic consumption rate

In the hyperbolic case, there is no exact analytical solution to equation (3.1). The general approach taken in the literature has been to use numerical methods to estimate the local partial pressure [11]. A drawback of this approach is the lack of a closed-form expression for the boundary of the anoxic radius, *r*_{n}. This differs from that in the constant model, as the reduced OCR in the hyperbolic formulation yields a longer diffusion distance. Denoting this modified anoxic radial distance by *a*(*p*)=0 at the anoxic edge *r*_{n} with a maximum value for *a*(*p*) at *r*_{o}. Close to the interface between oxygenated and hypoxic regions, the oxygen curve falls to zero with a zero gradient, which can be readily approximated by a quadratic parabola-like fit. If we define *a*_{o} as the maximum OCR, then a suitable approximation for *a*(*p*) is given by
*p*(*r*_{o})=*p*_{o}. From this,

### 3.2 Cylindrical symmetry

Cylindrically symmetric reaction–diffusion equations are often encountered in Krogh-type models to simulate oxygen flow from a blood vessel [8,9,13] under the reasonable assumption that vessels secrete oxygen perpendicular to their walls. In these situations, they present a reasonable approximation of a length of blood vessel running through a tumour. Such a situation is illustrated in figure 2. In this case, the steady-state oxygen distribution satisfies
*r*≤*r*_{c}, where *r*_{c} is the diffusion limit from the vessel.

#### 3.2.1 Constant consumption rate

With a constant consumption rate, *a*_{o}, an exact solution is possible. For a vessel of radius *r*_{o} with boundary conditions *p*(*r*_{c})=*p*′(*r*_{c})=0 and *p*(*r*_{o})=*p*_{o}, the solution to (3.8) is given by
*r*_{c} but this can be readily estimated using root-finding methods when equation (3.9) is equal to zero.

#### 3.2.2 Hyperbolic consumption rate

Analysis analogous to the spherical geometry can be applied in the hyperbolic case, with the caveat that the substitution term is modified, as in this case diffusion is out from the vessel and must be zero at the radius *p*(*r*_{o})=*p*_{o} yields an expression that is analytically tractable, but unwieldy. The first derivative of the solution is given by
*r* subject to *p*(*r*_{o})=*p*_{o} so that

### 3.3 Simulations and comparison with numerical solutions

In both geometries, the oxygen distributions resulting from both models were simulated with the maximum consumption rate parameter *a* varying over an order of magnitude so that 2.5×10^{−7}≤*a*≤1.5×10^{−6} m^{3} kg^{−1} s^{−1}, corresponding to the range 7.50–45.45 mm Hg *O*_{2} s^{−1}. For the spherical geometry, the external partial pressure was taken as 100 mm Hg. In the cylindrical geometry, both ‘arterial’ (100 mm Hg) and ‘venous’ (40 mm Hg) source partial pressures were considered and vessel radius was fixed at *r*_{o}=5 μm. These values were chosen as illustrations of the typical case, but in tumours chaotic vasculature often means significant variation in these values and can result in lower oxygen partial pressures [31], so for the cylindrical case a chaotic partial pressure of 20 mm Hg was also simulated. A value for the diffusion constant close to that of water (*D*=2×10^{−9} m^{2} s^{−1}) [8,9] was also used.

For each set of simulation parameters, the oxygen tension profile and diffusion length (anoxic radius in the spherical geometry, diffusion distance in the cylindrical geometry) were compared and the differences between both models for the given parameter set computed. The analytical approximations were contrasted with the numerical solutions to gauge how well they described the underlying behaviour of these models. As the value of *K*_{m} in the literature is typically ≤1 mm Hg [24,26], a value of 1 mm Hg was used in these simulations. To examine the sensitivity of the model to this parameter, simulations using values of 0 mm Hg≤*K*_{m}≤10 mm Hg were examined.

## 4. Results

### 4.1 Validity of the analytical approximation

To estimate the validity of our analytical approximation, we compared results from employing the approximation functions in both geometries with simulations using boundaries estimated from the approximation functions. Typical comparisons are shown in figure 3. In spherical geometry, the agreement is exceptionally good. Differences between the numerical simulation and analytical approximation were almost negligible, with only small disagreement at very low oxygen partial pressures (Δ*p*<0.24 mm Hg), indicating that the spherical approximation describes the solution to equation (3.1) well (all errors ≪1 mm Hg) and yields the behaviour expected from the equation. The agreement of the hyperbolic function with the approximation outlined in equation (3.5) is shown in figure 4.

In the cylindrical case, the analytical approximation systematically overestimated the oxygen partial pressure by a slight amount for each simulated case as can be seen in figure 3. This may be due to the two-term quadratic truncation introducing an error. The mean error was typically 0.3±1 mm Hg. Despite this systematic slight overestimation, the profile agreement is still high and the functional behaviour consistent. This indicates that the analytical approximation describes the function behaviour well in both geometries.

### 4.2 Comparison of constant and hyperbolic consumption rate models

Table 1 shows the estimated difference in boundary location between both models for a spherical geometry in the domain *r*_{n}, is of the order of tens of micrometres as shown in table 1. Despite the increased oxygen penetration and reduced anoxic radius in the hyperbolic model however, the oxygen profile is largely consistent between both models for spherical geometry. An example of this is illustrated in figure 5*a*, where the average difference between both profiles is only 0.35±0.60 mm Hg in the region *r*≤*r*_{o}) is considered, differences between both models are effectively negligible, being typically <0.5 mm Hg.

Similar trends are seen in the case of a cylindrical geometry, as shown in table 2. For the typical example illustrated in figure 5*b*, the average difference between the two models was 2.17±1.25 mm Hg over the entire profile, and a mean difference of 0.13±0.12 mm Hg in the region *r*_{c} is relatively small in this geometry at around 10 μm. However, over the entire region 0≤*r*≤*r*_{c}, differences between the constant consumption and hyperbolic consumption rate models are more pronounced in the cylindrical geometry.

Simulations of the hyperbolic model in both geometries were run with various values of *K*_{m} to investigate what effect this had on the expected oxygen profile. When *K*_{m}=0 mm Hg, the hyperbolic model reduces to the constant consumption case as expected. For higher values of *K*_{m}, the divergence between the constant and hyperbolic consumption cases increases with oxygen penetrating deeper and an increased overall oxygen profile. This behaviour is shown in figure 5 with *a*=5×10^{−7} m^{3} kg^{−1} s^{−1}. At *K*_{m}=10 mm Hg, an order of magnitude above literature values, oxygen diffuses approximately 55% further through the spheroid than in the constant consumption case. In the cylindrical case at *K*_{m}=10 mm Hg, the diffusion distance is 30% greater than in the constant consumption case. This significantly changes the oxygen profile in both geometries as well as the diffusion distance with markedly different values expected for different values of *K*_{m}. However, as literature values to date indicate that *K*_{m}<1 mm Hg, this would suggest that the expected oxygen profile in all geometries should be very close to that of the constant consumption rate case. Consequently, it is unlikely that this issue would ever be a significant factor when oxygen is the substrate but may be an issue for applying these dynamics to substances other than oxygen.

## 5. Discussion

While the differences in oxygen profile between the constant and hyperbolic consumption rate models are relatively small, the increased oxygen diffusion distance in the hyperbolic case may have an impact when estimating the extent of the hypoxic region. In the regions of difference between the constant and hyperbolic rate models (

In spherical geometry, the divergence between the constant and hyperbolic consumption model is generally relatively minor, as can be seen in table 1. There are however situations where this is not the case; for example, spheroids with low consumption rates and a radius of less than 500 *μ*m, such as the first entry in table 1 for a spheroid with radius 500 μm and OCR of 2.5×10^{−7} m^{3} kg^{−1} s^{−1}. In this case, the difference between *r*_{n} and *K*_{m}, it might be possible to use stained spheroids to investigate this further. One option would be spheroids stained with EF5, which binds to oxygenated regions between 0.8 and 10 mm Hg [32]. This would manifest in small differences between the expected location of the stained boundary in both models, maximal at *p*≃0.8 mm Hg. In both spherical and cylindrical geometries, the projected differences are typically *a* is of the same order of magnitude as the uncertainties in estimating the stained boundary [10], rendering direct comparison extremely difficult. Conversely, this implies that at current experimental resolution, the differences between the two models are negligible if maximal value of *a* can be accurately deduced. This suggests that the constant OCR model is a suitable model for this geometry regardless of whether the OCR is governed by hyperbolic kinetics or not.

In cylindrical geometry, differences between the constant and hyperbolic consumption rate models were small but still more pronounced than in the spherical case. It is important to note that differences between the two models in cylindrical geometry may be considerably less than that presented here as the analytical approximation systematically overestimates the oxygen profile. Hence, it follows that differences between the constant and hyperbolic consumption rate models are even less than the values reported here. More experimental data are required to investigate this further. The models derived in this work for the cylindrical geometry are Krogh-type models, which have been commonly used to model vessel oxygenation and produce results in good agreement with clinical data [6,8,9]. This approach is not by any means the sole method of estimating oxygen diffusion, and other authors have used numerical methods, including Green’s function methods [14] to solve the reaction–diffusion equations for vessel geometry. A drawback of this method is the difficulty in assigning appropriate boundary conditions [29], a difficulty which appears to some extent in all numerical approaches to oxygen modelling. One advantage of the models proposed in this work is that the boundaries for oxygen diffusion can be readily estimated, circumventing this issue and allowing the establishment of explicit relationships between the OCR and the resultant oxygen diffusion distance. In this work, the illustrative examples show the simplest case of vessels with a constant value for *p*_{o} along the length. One limitation of the cylindrical model is that is does not take blood flow or number of erythrocytes in the capillary into account; this could have a significant effect, as the oxygen partial pressure may decrease markedly along a vessel length. If the variation rate of oxygen pressure along the vessel length can be estimated with the relevant perfusion parameters, then it should be possible to extend the simple model presented here to consider this factor.

To gauge sensitivity to the parameter *K*_{m}, simulations were run from zero to an order of magnitude above literature values (0≤*K*_{m}≤10 mm Hg). Values of *K*_{m} far above literature values greatly increased the diffusion differences between both geometries, and increased the overall oxygen profile significantly. As *K*_{m} tended to zero, hyperbolic consumption rate models in both geometries reduced to constant consumption models. As literature values are typically *K*_{m}≤1 mm Hg for oxygen, the effects of higher values can likely be neglected.

Whether constant OCR or the hyperbolic approximation models are used, the effect of *a*_{o} on oxygen diffusion distance is important and worthy of discussion. Intuitively, we might expect that higher values of maximum consumption rate *a*_{o} mean that oxygen is depleted rapidly and as a consequence the oxygen diffusion distance is diminished. The simulations in this work reiterate this in both geometries; in the spherical configuration, increasing *a*_{o} corresponded to decreasing values for the viable spherical rim *r*_{c} and increased values for the hypoxic radius *r*_{n}. In cylindrical vessel-like geometry, increasing *a*_{o} acted to reduce the diffusion distance *r*_{l}, as illustrated in figure 6. Unlike many numerical approaches, the relationship between *a*_{o} and diffusion distance is explicitly derived in this work. This analysis suggests that OCR has a marked effect on tumour hypoxia in both vascular and avascular situations—for the former, high values of *a*_{o} would act to reduce the distance oxygen can travel in the tissue, potentially increasing the prevalence of hypoxic regions. In the avascular case, high values of *a*_{o} would result in increased central anoxia. While beyond the scope of this work, this aspect is worthy of further investigation.

The results of this analysis would suggest that in general, constant consumption rate type models yield a broadly similar oxygen profile to hyperbolic-type models and are relatively easy to implement, even for large datasets. There are however situations where significant differences can develop between both approaches, and the forms derived here allow rapid quantification of the resultant oxygen profiles and diffusion distances for either form of reaction term used. Although broadly similar for most realistic parameter values, there are some inherent differences between constant oxygen consumption models and hyperbolic consumption forms, chiefly that the hyperbolic case always results in a deeper diffusion depth than for a constant consumption model with the same parameters. This difference is intuitively grasped by the consideration that OCR falls with decreasing oxygen tension in such models and as a result oxygen can penetrate further at low oxygen levels than in a comparable constant consumption case. This analysis suggests that such differences are usually minor and not typically resolvable through clinical means, but also that there are situations where a measurable difference between both approaches could arise. In the formulation established in this study, the explicit relationship between maximal OCR *a* and the diffusion distances in all cases can be estimated from first principles, allowing rapid quantification of oxygen boundaries in all cases.

## 6. Conclusion

This analysis suggests that since the effective *K*_{m} constant for oxygen consumption is so small, only minor differences in diffusion limit (≃10 μm) and local oxygen partial pressure (Δ*p*≤1 mm Hg in spherical geometry, Δ*p*≃2 mm Hg in cylindrical geometry) would be expected between the constant and hyperbolic consumption rate models for the typical range of literature values of *K*_{m}. There are also cases where differences between both forms can be significant, and the forms outlined in this work allow rapid first principles estimate of the respective oxygen profiles and diffusion distance for both constant consumption rate and hyperbolic consumption rate model forms in spherical and cylindrical geometry.

In general, the simpler constant consumption form can be used in both spherical and cylindrical geometries for most situations. The models outlined also allow quick identification of situations where both models are expected to produce significantly different results, and suggestions even in these limited situations the oxygen profiles for both approaches are broadly similar. In oxygen modelling situations where hyperbolic kinetics are required, the analytical steady-state approximations outlined in this work could be used readily, allowing rapid estimation of the oxygen diffusion limit and projected oxygen profile.

## Data Accessibility

The full solution of equation (3.12) can be found in Mathematica Format uploaded as supplementary material.

## Funding statement

D.R.G. and M.P. are funded by Cancer Research UK and the Medical Research Council. A.G.F. is supported by the EPSRC and Microsoft Research, Cambridge through grant no. EP/I017909/1 (www.2020science.net).

- Received June 3, 2014.
- Accepted July 30, 2014.

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.