A Mathematical Dissection of the Adaptation of Cell Populations to Fluctuating Oxygen Levels

Abstract

The disordered network of blood vessels that arises from tumour angiogenesis results in variations in the delivery of oxygen into the tumour tissue. This brings about regions of chronic hypoxia (i.e. sustained low oxygen levels) and regions with alternating periods of low and relatively higher oxygen levels, and makes it necessary for cancer cells to adapt to fluctuating environmental conditions. We use a phenotype-structured model to dissect the evolutionary dynamics of cell populations exposed to fluctuating oxygen levels. In this model, the phenotypic state of every cell is described by a continuous variable that provides a simple representation of its metabolic phenotype, ranging from fully oxidative to fully glycolytic, and cells are grouped into two competing populations that undergo heritable, spontaneous phenotypic variations at different rates. Model simulations indicate that, depending on the rate at which oxygen is consumed by the cells, dynamic nonlinear interactions between cells and oxygen can stimulate chronic hypoxia and cycling hypoxia. Moreover, the model supports the idea that under chronic-hypoxic conditions lower rates of phenotypic variation lead to a competitive advantage, whereas higher rates of phenotypic variation can confer a competitive advantage under cycling-hypoxic conditions. In the latter case, the numerical results obtained show that bet-hedging evolutionary strategies, whereby cells switch between oxidative and glycolytic phenotypes, can spontaneously emerge. We explain how these results can shed light on the evolutionary process that may underpin the emergence of phenotypic heterogeneity in vascularised tumours.

Analysis of Evolutionary Dynamics for a Simplified Model

In order to obtain a detailed analytical description of the evolutionary dynamics of the two cell populations, we can consider a simplified scenario whereby the ODE for S(t) is decoupled from the system of non-local parabolic PDEs (4). In particular, we let the evolution of the oxygen concentration S(t) be governed by the following Cauchy problem

$$\begin{aligned} \left\{ \begin{array}{ll} {\frac{\mathrm{d}S}{\mathrm{d}t} = I(t)-\Theta S},\\ S(0) = S^0 \ge 0, \end{array} \right. \quad t \in (0,\infty ), \end{aligned}$$

(17)

where the effects of oxygen consumption and oxygen decay are both encapsulated in the parameter $\Theta >0$. Moreover, to facilitate analysis, we extend the interval [0, 1] to $\mathbb {R}$ and re-define the population-level quantities accordingly, i.e. we use the definitions

$$\begin{aligned} \rho _{H}(t) := \int _{\mathbb {R}} n_{H}(x,t) \; \mathrm{d}x, \quad \rho _{L}(t) := \int _{\mathbb {R}} n_{L}(x,t) \; \mathrm{d}x ,\quad \rho (t) := \rho _H(t) + \rho _L(t) \end{aligned}$$

and

$$\begin{aligned} \mu _{i}(t) := \frac{1}{\rho _i(t)} \int _{\mathbb {R}} x \, n_{i}(x,t) \; \mathrm{d}x, \quad \sigma ^2_{i}(t) := \frac{1}{\rho _i(t)} \int _{\mathbb {R}} x^2 \; n_i(x,t) \; \mathrm{d}x - \mu _i^2(t) \end{aligned}$$

with $i \in \{H,L\}$. Finally, in agreement with much of the previous work on the mathematical analysis of the evolutionary dynamics of continuous traits, which relies on the simplifying assumption that population densities are Gaussians (Rice 2004), we consider initial conditions of the form

$$\begin{aligned} n_i(x,0) = \rho _i^0 \sqrt{\frac{v_i^0}{2\pi }} \exp \left[ -\frac{v_i^0}{2} (x-\mu _i^0)^2 \right] , \; \text { where } \; \rho _i^0, v_i^0\in \mathbb {R}_{>0} \; \text { and } \; \mu _i^0\in \mathbb {R}. \end{aligned}$$

(18)

This allows us to use the result established by Proposition 1, which can be proved through the method that we previously employed in Ardaševa et al. (2020).

Proposition 1

Under assumptions (6) and (7), the system of non-local PDEs (4) posed on $\mathbb {R} \times (0,\infty )$ and subject to the initial condition (18) admits the exact solution

$$\begin{aligned} n_i(x,t)= \rho _i(t) \sqrt{\frac{v_i(t)}{2 \pi }} \exp \left[ -\frac{v_i(t)}{2} (x-\mu _i(t))^2 \right] \quad \text {for} \quad i \in \left\{ H, L \right\} , \end{aligned}$$

(19)

with the population size, $\rho _i(t)$, the mean phenotype, $\mu _i(t)$, and the inverse of the phenotypic variance, $v_i(t) = 1/\sigma ^{2}_i(t)$, being solutions of the Cauchy problem

$$\begin{aligned} \left\{ \begin{array}{ll} \displaystyle {\frac{\mathrm{d}v_i}{\mathrm{d}t} = 2 \left( h(S) - \beta _i v_i^2\right) },\\ \displaystyle {\frac{\mathrm{d}\mu _i}{\mathrm{d}t} = \frac{2h(S)}{v_i} \left( \varphi (S) - \mu _i\right) },\\ \displaystyle {\frac{\mathrm{d}\rho _i}{\mathrm{d}t} = \left( F_i(S,v_i,\mu _i) - d \rho \right) \rho _i},\\ v_i(0) = v^0_i, \quad \mu _i(0) = \mu ^0_i, \quad \rho _i(0) = \rho ^0_i,\\ \rho := \rho _H + \rho _L, \end{array} \right. \qquad \text {for } \; i \in \left\{ H, L \right\} , \end{aligned}$$

(20)

where

$$\begin{aligned} F_i(S,v_i,\mu _i) := \gamma \, g(S) - \frac{h(S)}{v_i} - h(S) \left( \mu _i - \varphi (S) \right) ^2. \end{aligned}$$

(21)

In the case where the inflow of oxygen is constant, i.e. the source term I(t) in the ODE (17) satisfies assumption (14), our main results are summarised by Theorem 1, where the functions g, $\varphi$ and h are defined according to (9) and (10), and we use the definitions

$$\begin{aligned} S^{\infty } := \frac{I_S}{\Theta }, \quad \rho ^{\infty }_L := \frac{ \gamma \, g(S^{\infty }) - \sqrt{h(S^{\infty }) \, \beta _L}}{d}, \quad \mu ^{\infty }_L := \varphi (S^{\infty }). \end{aligned}$$

(22)

Theorem 1

Under assumptions (5)–(11) and the additional assumption (14), the solution of the system of non-local PDEs (4) posed on $\mathbb {R} \times (0,\infty )$, subject to the initial condition (18) and complemented with the Cauchy problem (17) is of the Gaussian form (19) and satisfies the following:

• (i) if

  $$\begin{aligned} \sqrt{h(S^{\infty }) \, \beta _L} \ge \gamma \, g(S^{\infty }) \end{aligned}$$

  then

  $$\begin{aligned} \lim _{t\rightarrow \infty }\rho _H(t)=0 \quad \text {and} \quad \lim _{t\rightarrow \infty }\rho _L(t)=0; \end{aligned}$$

• (ii) if

  $$\begin{aligned} \sqrt{h(S^{\infty }) \, \beta _L} < \gamma \, g(S^{\infty }) \end{aligned}$$

  then

  $$\begin{aligned} \lim _{t\rightarrow \infty }\rho _H(t)=0, \quad \lim _{t\rightarrow \infty }\rho _L(t)= \rho ^{\infty }_L \end{aligned}$$

  and

  $$\begin{aligned} \lim _{t\rightarrow \infty }\mu _L(t) = \mu ^{\infty }_L, \quad \lim _{t\rightarrow \infty } \sigma ^2_L(t) = \sqrt{\frac{\beta _L}{h(S^{\infty })}}. \end{aligned}$$

In the case where the inflow of oxygen undergoes periodic oscillations, i.e. the source term I(t) in the ODE (17) satisfies assumption (15) along with the additional assumption

$$\begin{aligned} I \in \mathrm{Lip}([0,\infty )), \end{aligned}$$

(23)

our main results are summarised by Theorem 2, where $\tilde{S}(t)$ is the unique non-negative T-periodic solution of the problem

$$\begin{aligned} \left\{ \begin{array}{ll} \displaystyle {\frac{\mathrm{d}\tilde{S}}{\mathrm{d}t} = I(t)-\Theta \tilde{S}, \quad t \in (0,T),}\\ \tilde{S}(0) = \tilde{S}(T), \end{array} \right. \end{aligned}$$

(24)

$\tilde{v}_i(t)$ is the unique real positive T-periodic solution of the problem

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle {\frac{\mathrm{d}\tilde{v}_i}{\mathrm{d}t} = 2 \left( h(\tilde{S}) - \beta _i \tilde{v}_i^2\right) , \quad t \in (0,T),}\\ \tilde{v}_i(0) = \tilde{v}_i(T), \end{array}\right. } \end{aligned}$$

(25)

$\tilde{\mu }_i(t)$ is the unique real T-periodic solution of the problem

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle {\frac{\mathrm{d}\tilde{\mu }_i}{\mathrm{d}t} = \frac{2 h(\tilde{S})}{\tilde{v}_i} \left( \varphi (\tilde{S}) - \tilde{\mu }_i\right) , \quad t \in (0,T),} \\ \tilde{\mu }_i(0) = \tilde{\mu }_i(T), \end{array}\right. } \end{aligned}$$

(26)

$\tilde{\rho }_i(t)$ is the unique real non-negative T-periodic solution of the problem

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle {\frac{\mathrm{d}\tilde{\rho }_i}{\mathrm{d}t} = \left( F_i(\tilde{S}, \tilde{v}_i,\tilde{\mu }_i) - d \tilde{\rho }_i \right) \tilde{\rho }_i, \quad t \in (0,T)}, \\ \\ \tilde{\rho }_i(0) = \tilde{\rho }_i(T), \end{array}\right. } \end{aligned}$$

(27)

and

$$\begin{aligned} \Lambda _i := \frac{1}{T} \int _{0}^{T} \frac{h(\tilde{S}(z))}{\tilde{v}_i(z)} \, \mathrm{d}z + \frac{1}{T} \int _{0}^{T} \left( \tilde{\mu }_i(z) - \varphi (\tilde{S}(z)) \right) ^2 h(\tilde{S}(z)) \, \mathrm{d}z \quad \text {for} \quad i \in \left\{ H,L\right\} . \end{aligned}$$

(28)

In (25)–(28), the functions g, $\varphi$ and h are defined according to (9) and (10). Moreover, the function $F_i$ in (27) is defined according to (21).

Theorem 2

Under assumptions (5)–(11) and the additional assumptions (15) and (23), the solution of the system of non-local PDEs (4) posed on $\mathbb {R} \times (0,\infty )$, subject to the initial condition (18) and complemented with the Cauchy problem (17) is of the Gaussian form (19) and satisfies the following:

• (i) if

  $$\begin{aligned} \min \left\{ \Lambda _H, \Lambda _L\right\} \; \ge \; \frac{\gamma }{T} \int _{0}^{T} g(\tilde{S}(t)) \, \mathrm{d}t \end{aligned}$$

  then

  $$\begin{aligned} \lim _{t\rightarrow \infty }\rho _H(t)=0 \quad \text {and} \quad \lim _{t\rightarrow \infty }\rho _L(t)=0; \end{aligned}$$

• (ii) if

  $$\begin{aligned} \min \left\{ \Lambda _H, \Lambda _L\right\} < \; \frac{\gamma }{T} \int _{0}^{T} g(\tilde{S}(t)) \, \mathrm{d}t \end{aligned}$$

  and

  $$\begin{aligned} i = \underset{k \in \left\{ H,L \right\} }{{\text {arg}}\,{\text {min}}}\; \Lambda _k, \quad j = \underset{k \in \left\{ H,L \right\} }{{\text {arg}}\,{\text {max}}}\; \Lambda _k, \end{aligned}$$

  then

  $$\begin{aligned} \rho _i(t) \rightarrow \tilde{\rho }_i(t), \quad \rho _j(t) \rightarrow 0 \quad \text {as } t \rightarrow \infty , \end{aligned}$$

  and

  $$\begin{aligned} \mu _i(t) \rightarrow \tilde{\mu }_i(t), \quad \sigma ^2_i(t) \rightarrow \frac{1}{\tilde{v}_i(t)} \quad \text {as } t \rightarrow \infty . \end{aligned}$$

Theorems 1 and 2 can be proved through methods similar to those that we employed in Ardaševa et al. (2020) and, therefore, their proofs are omitted here.