Oxygen consumption dynamics in steady-state tumour models

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

νmax2.1

where s is the substrate concentration, ${\nu }_{max}$ is the maximum reaction rate and K_m is substrate level where ν is half ${\nu }_{max}$.

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

$$0=\frac{D}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2\frac{\mathrm{d}p}{\mathrm{d}r}\right)-\mathit{\Omega }a(p)$$3.1

in the region r_n<r<r_o, where Ω=3.0318×10⁷ mm Hg kg m⁻³ 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 prioribut 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

$$p(r)=p_o+\frac{a_o\mathit{\Omega }}{6D}(r^2-r_o^2+2r_n^3(\frac{1}{r}-\frac{1}{r_o}))=\frac{a_o\mathit{\Omega }}{6D}(r^2+\frac{2r_n^3}{r}-3r_n^2)$$3.2

for 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_l=\sqrt{\frac{6D{p}_o}{\mathit{\Omega }{a}_o}}.$$3.3

It is important to note that 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=r_o(\frac{1}{2}-\cos \left(\frac{\arccos (1-2r_l^2/r_o^2)-2\pi }{3}\right)).$$3.4

For spheroids sufficiently small to have no central anoxia, 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 $r_n^{\ast }$, we can establish a phenomenological approximation which captures the behaviour of the function and satisfies the conditions that 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

$$a(p)=\frac{a_{o}p(r)}{p(r)+K_m}\approx \frac{a_o{(r-r_n^{\ast })}^2}{{(r-r_n^{\ast })}^2+k_r},$$3.5

where

$$k_r={(r_o-r_n^{\ast })}^2\frac{K_m}{p_o}.$$3.6

Substituting (3.5) into (3.1) with the same boundary conditions as the constant consumption rate case, $p(r_n^{\ast })=0$ and $p^{\prime }(r_n^{\ast })=0$, we obtain the analytic solution

$$\begin{array}{rl}p(r)& =\frac{a_o\mathit{\Omega }}{6D}(r^2-3{r_n^{\ast }}^2+\frac{2{r_n^{\ast }}^3-6k_r{r}_n^{\ast }}{r}+6k_r-3\sqrt{k_r}\pi r_n^{\ast }+6\sqrt{k_r}{r}_n^{\ast }\arctan \left(\frac{\sqrt{k_r}}{r-r_n^{\ast }}\right)\\ & +\frac{(6\sqrt{k_r}{r_n^{\ast }}^2-6k_r^{3/2})\arctan ((r-r_n^{\ast })/\sqrt{k_r})}{r}+(\frac{6k_r{r}_n^{\ast }}{r}-3k_r)\ln \left(\frac{{(r-r_n^{\ast })}^2+k_r}{k_r}\right)).\end{array}$$3.7

To estimate $r_n^{\ast }$ in this case, (3.8) is solved subject to p(r_o)=p_o. From this, $r_n^{\ast }$ can be found to any required degree of precision from the resulting expression, allowing determination of this boundary. It is important to note that oxygen can diffuse further in the hyperbolic formulation, and subsequently the anoxic region is decreased relative to the constant consumption case so that $r_n^{\ast }\le r_n$. In the case of spheroids sufficiently small to have no central anoxia, $r_n^{\ast }=0$.

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

$$0=\frac{D}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left(r\frac{\mathrm{d}p}{\mathrm{d}r}\right)-a(p)\mathit{\Omega }$$3.8

for all 0≤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

$$p(r)=p_o+\frac{a_o\mathit{\Omega }}{4D}(r^2-r_o^2-2r_c^2\ln \frac{r}{r_o}).$$3.9

Unlike the spherically symmetric case, there is no corresponding analytical expression for the diffusion distance 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 $r_c^{\ast }$. In this case, oxygen should diffuse further from the vessel than in the constant consumption rate case. This approximation is given by

$$a(p)=\frac{a_{o}p(r)}{p(r)+K_m}\simeq \frac{a_o{(r_c^{\ast }-r)}^2}{{(r_c^{\ast }-r)}^2+l_r},$$3.10

where

$$l_r={(r_c^{\ast }-r_o)}^2\frac{K_m}{p_o}.$$3.11

Solving (3.8) subject to the boundary conditions $p(r_c^{\ast })=p^{\prime }(r_c^{\ast })=0$ and p(r_o)=p_o yields an expression that is analytically tractable, but unwieldy. The first derivative of the solution is given by

$$p^{\prime }(r)=\frac{a_o\mathit{\Omega }}{2D}(r+\frac{2r_c^{\ast }\sqrt{l_r}\arctan ((r_c^{\ast }-r)/\sqrt{l_r})-l_r\ln ({(r_c^{\ast }-r)}^2/l_r)-{r_c^{\ast }}^2}{r}).$$3.12

To find $r_c^{\ast }$, the modified diffusion distance, equation (3.12) is integrated with respect to r subject to $p(r_c^{\ast })=0$. Equation (3.12) can be solved analytically, producing a long expression which is not shown here for the sake of brevity. The full solution and code is supplied in the electronic supplementary material for reference. The resulting form is then evaluated with condition p(r_o)=p_o so that $r_c^{\ast }$ may be found to any desired degree of accuracy.

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⁻⁷≤a≤1.5×10⁻⁶ m³ kg⁻¹ s⁻¹, corresponding to the range 7.50–45.45 mm Hg O₂ s⁻¹. 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⁻⁹ m² s⁻¹) [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.