Environmental stress level to model tumor cell growth and survival

1 Introduction

As cancerous diseases still cause many deaths in human population, there is a high need of improving the understanding of tumor growth dynamics and treatment strategies. For that purpose and for treatment optimization, mathematical approaches can help with adequate quantitative descriptions and predictions. Numerous mathematical models and methods are already available (see [1, 2]), including differential equations of all types, stochastic models, phenomenological models as well as more mechanistic ones. Despite good new experimental techniques, e.g. measuring tumor growth also in vivo [3, 4, 5], it is still difficult to quantify growth by determining parameter values in realistic situations without harming a patient or affecting a potential treatment. To get a better quantitative understanding of the underlying processes, investigating tumor cell cultures and even single cells under well-controlled conditions provides a promising way to get at least parts of the necessary information (see e.g. [6]). But even under laboratory conditions, parameters can underlie variations and cells may behave individually to some extent. To consider this, uncertainty quantification can be an adequate tool. Having such in vitro

experiments, it is much easier to separate and measure different environmental influences on the tumor growth. The probably most relevant environmental influences are the availability of nutrients and oxygen. For modeling purposes, one often focuses on the limiting nutrient, without specifying it in detail. Nutrient deprivation causes “stress” to any cells, but other factors may cause stress to the cells, as well. Thus, it is worth to consider the concept of stress as a relevant factor for modeling growth and death of tumor cells.

To approach this problem, we suggest the introduction of an “environmental stress level” as an immeasurable auxiliary variable, which quantifies to what extent viable cells would get stressed out under certain environmental conditions. These may include quantities like nutrient saturation, oxygen level, drug concentrations or mechanical forces. On the one hand, this environmental stress level can e.g. inhibit cell growth and promote death by inducing necrosis for present cells. On the other hand, it is assumed to have an influence on cell movement, as cells favor areas with a low environmental stress level, which may cause metastasis formation, as well. In this work, this approach is presented in a simplified spatially homogeneous setting, where only the nutrient saturation effects the environmental stress level. Nevertheless, the presented model can be easily extended to more environmental variables, which provides several advantages from a modeling point of view. Even in this very simplified setting, model comparison shows a slight preference for the newly introduced model.

For comparison, we present two systems of ordinary differential equations, which model tumor cell population dynamics respectively with and without considering an environmental stress level. To do so, time-resolved in vitro data (from [7]

) are used to perform model parameter calibration with Bayesian inversion. The latter approach provides a quantitative estimate of the uncertainty in the estimated parameters and hence brings more information than a deterministic inversion method, at the price of a higher computational cost. In the Bayesian inversion setting several methods are available for sampling. Among the most popular ones are Markov Chain Monte Carlo (MCMC) methods

[8, Ch. 6-7], which are very simple to implement, but they need ad hoc tuning when the dimension of the parameter space is high or for complex posterior distributions [9]. In the last decade, particle based methods have gained popularity for their higher efficiency and robustness compared to simple MCMC, especially in the presence of high dimensional parameter spaces and multimodal posterior distributions [10]. Most of these algorithms are population based, which means they deal with a collection of samples in every iteration. In this category are sequential importance sampling and resampling (annealed importance sampling [11] and Sequential Monte Carlo [12, 13]) and population MCMC [14, 15]

. To calibrate the unknown parameters of the presented models, we use Sequential Monte Carlo (SMC). It is worth mentioning that the ensemble Kalman filter

[16] can also be used to solve the inverse problem. It can be viewed as an approximate SMC with Gaussian approximation of the posterior [17]

. However, the rigorous analysis of this method is confined for the moment to linear forward operators

[18, 17]

. Eventually, the calibration results are used to compare the models via a validation metric and the Bayes factor. For this purpose, the chosen SMC method is advantageous opposed to other algorithms (e.g. MCMC), since its structure provides easy access to the Bayes factor

[10].

The article is structured as follows: The upcoming Section 2 explains the investigated mathematical models with the underlying experimental setup and provides background about the applied methods for uncertainty quantification, parameter calibration and model comparison. The results from the model calibrations are summarized and discussed in Section 3. Finally, Section 4 concludes the presented results and the use of the novel modeling approach. At the end of this article, more details about mathematical analysis of the models and the calibration results are enclosed in a supplementary section.

2 Materials and method

Section 2.1 presents different models for population dynamics of tumor cells and introduces the idea of using a stress level to model environmental influences on the cells collectively. To investigate the effect of nutrient changes on viable cells, we consider a special case of the models. The underlying biological setting and the measurement methods for the corresponding in vitro experiments are explained in Section 2.2. With these data we calibrate the unknown parameters of the models using Bayesian inversion and Sequential Monte Carlo methods. Section 2.3 provides the theoretical background behind the algorithms. Finally, we present the validation metric and the Bayes factor in Section 2.4, which are used to investigate and compare the resulting model calibrations.

2.1 Mathematical modeling

The presented ODE models describe the basic growth and death dynamics of tumor cells in a spatially homogeneous, avascular environment under consideration of the present nutrient saturation. The time-dependent variables are the density of viable tumor cells $V = V (t)$ , the nutrient saturation $S = S (t)$ , and the environmental stress level $η = η (t)$ at day $t \geq 0$ .

We consider two possibilities to model the influence of the nutrient saturation on the cells: Model $M_{S}$ assumes that the nutrient saturation directly affects the growth/death rates of the cells, while in model $M_{η}$ the nutrient supply influences the cell population indirectly by changing the environmental stress level, which itself has an effect on the viability of the cells. The nutrient saturation is bounded and normalized: $S (t) \in [0, 1]$ . The bounds represent the complete absence of nutrients ( $S = 0$ ) and a nutrient supply, which creates optimal growth conditions for the cells ( $S = 1$ ). In the latter case, both models can be reduced to a simplified model, which is independent of the nutrient saturation $S$ . We denote this model by $M_{opt}$ (“opt” short for “optimal” growth conditions).

A detailed description of each model follows in the consecutive paragraphs. All models are biologically reasonable in a first check, as they have been mathematically analyzed in terms of positivity and boundedness of the solutions as well as steady states and their stability. The last paragraph of this section provides a summary of the variables, parameters, and mathematical features of the models.

What if the cells find optimal nutrient conditions?

External nutrients serve as an energy source for cell proliferation. Growth can happen with a maximal possible proliferation rate $β$ , which is virtually reached under optimal nutrient conditions ( $S = 1$ ). Under these circumstances, cell death is assumed to occur only for nutrient-independent reasons (e.g. dying of old age) with a constant rate $λ$ . These dynamics can be modeled with a combination of a generalized logistic growth term (first proposed in [19]) and an exponential death term. The population growth is limited by the carrying capacity $K$ of the biological system, taking into account e.g. limited space but not nutrient shortage. In particular, cell growth is inhibited in a cell density dependent manner, which is known as “proliferation contact inhibition” [20]. The strength of this phenomenon is specific for the cell type and is modeled by the shape parameter $m$ of the logistic growth. A small value for $m$ indicates strong contact inhibition, i.e. reduction of proliferation already starts at small cell densities. Therefore, $1 / m$ can be interpreted as the strength of contact inhibition. Overall, this results in the following initial value problem:

S = 1 \Rightarrow ⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} ˙ V & = β V (1 - {(\frac{V}{K})}^{m}) - λ V, V (0) & = V_{0} . \end{matrix}

(

M_{opt}

)

As a feature of the pure logistic term, population growth for very small populations (i.e. $V ≪ K$ ) can be approximated by exponential growth with rate $β$ . It is a reasonable assumption that, independently of the initial population size $V_{0}$ , the cell population’s size actually increases over time under optimal nutrient conditions, which translates to a parameter relation: $β > λ$ . With this constraint, the ODE can be rewritten as a purely logistic growth term with a new rate $~ β$ and capacity $~ K$ :

	$˙ V$

Therefore, under optimal nutrient conditions, the cell population actually grows logistically with

	$``net'' growth rate: ~ β$	$= (β - λ) \in (0, β)$
	$and ``net'' carrying capacity: ~ K$	$= K {(1 - \frac{λ}{β})}^{\frac{1}{m}} = K {(\frac{}{β - λ} β)}^{\frac{1}{m}} \in (0, K) .$

For $t \geq 0$ the known analytical solution of such an initial value problem, i.e. for model $M_{opt}$ , is

	$V (t)$	$= \frac{V_{0} ~ K}{{(V_{0}^{m} + ({~ K}^{m} - V_{0}^{m}) \cdot e^{- ~ β m t})}^{\frac{1}{m}}}$
		$= {(\frac{(V_{0} K)^{m} (β - λ)}{β V_{0}^{m} + ((β - λ) K^{m} - β V_{0}^{m}) \cdot e^{- (β - λ) m t}})}^{\frac{1}{m}} .$		(2.1)

What if suboptimal nutrient conditions directly affect the cell dynamics?

Exposing the cells to suboptimal nutrient conditions ( $0 \leq S < 1$ ) can affect the cells’ metabolism by decelerating proliferation or even lead to cell death by starvation. To include this, an additional death term is needed, whose rate as well as the proliferation rate are scaled appropriately to the nutrient supply.

A time-dependent nutrient supply can be described by adding an additional differential equation to the system: $˙ S = g (S, V, t)$ with $S (0) = S_{0}$ . The reaction term $g (S, V, t)$ can e.g. include terms for nutrient consumption by viable cells or an external nutrient source. In the considered in vitro experiments, whose measurements are used to calibrate the model later on, the nutrient concentration is maintained approximately constant. Therefore, all presented models assume

˙ S = g (S, V, t) = 0 \Rightarrow S (t) = const. = S_{0} \in [0, 1] \forall t \geq 0 .

The present nutrient saturation $S_{0}$ contributes to the cells’ reproductivity as well as to their ability to survive. If there are no nutrients available, the cells die from starvation with a maximal possible rate $λ_{st}$ , assumed to happen in an exponential manner. With increasing nutrient supply, the starvation rate $λ_{st}$ decreases and the proliferation rate $β$ is up-regulated. Adapting the previous model $M_{opt}$ accordingly yields the more general initial value problem

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ \begin{matrix} ˙ V & = δ^{+} (S_{0}) \cdot β      = β_{S} V (1 - {(\frac{V}{K})}^{m}) - (λ + δ^{-} (S_{0}) \cdot λ_{st})      = λ_{S} V, V (0) & = V_{0}, \end{matrix}

(

M_{S}

)

with appropriate nutrient-dependent functions $δ^{+} (S)$ and $δ^{-} (S)$ , scaling the total growth and death rates:

β_{S} = δ^{+} (S_{0}) \cdot β and λ_{S} = (λ + δ^{-} (S_{0}) \cdot λ_{st}) .

The scaling functions $δ^{\pm} (S)$ need to fulfill certain criteria to be biologically meaningful in our model setting. First, the cells should not die from starvation but proliferate with the maximal possible rate under optimal nutrient conditions. In contrast, without nutrients the cells should not be able to reproduce but starve with a maximal possible rate. These features translate to the properties

δ^{+} (S) \cdot β ⎧ ⎨ ⎩ \begin{matrix} \to β & for S \to 1, \in (0, β) & for 0 < S < 1, \to 0 & for S \to 0, \end{matrix} and δ^{-} (S) \cdot λ_{st} ⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} \to 0 & for S \to 1, \in (0, λ_{st}) & for 0 < S < 1, \to λ_{st} & for S \to 0 . \end{matrix}

Second, increasing the nutrient saturation should promote growth and reduce starvation: $δ^{+} (S_{0})$ resp. $δ^{-} (S_{0})$ need to be monotonically increasing resp. decreasing functions bounded by $[0, 1]$ . Therefore, a possible choice are Hill type functions with Hill coefficient $k = 2$ , i.e.

δ^{+} : S \mapsto \frac{S^{2}}{S_{thr}^{2} + S^{2}} and δ^{-} : S \mapsto 1 - δ^{+} (S) = 1 - \frac{S^{2}}{S_{thr}^{2} + S^{2}},

(2.2)

where $0 < S_{thr} ≪ 1$ denotes the nutrient threshold, for which the cells proliferate/starve with half-maximal rate. For simplicity, we do not assume any hysteresis effect, motivating the relation $δ^{-} (S) = 1 - δ^{+} (S)$ . The scaling functions $δ^{\pm} (S)$ can be interpreted as influence functions, which describe how tolerant viable cells are to nutrient changes. In this context, the parameter $S_{thr}$ can be seen as a nutrient sensitivity threshold. Figure 1 shows a qualitative plot of the behavior of the scaling functions.

Figure 1: Qualitative behavior of the Hill type scaling functions $δ^{+} (S)$ and $δ^{-} (S)$

It can be shown (for calculations see Section 5.1 in the supplement) that for $t \geq 0$ a well-defined analytical solution of problem ( $M_{S}$ ) is

V (t) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} V_{0} K \cdot {⎛ ⎜ ⎝ \frac{β_{S} - λ_{S}}{β_{S} V_{0}^{m} + ((β_{S} - λ_{S}) K^{m} - β_{S} V_{0}^{m}) e^{- (β_{S} - λ_{S}) m t}} ⎞ ⎟ ⎠}^{\frac{1}{m}} & if β_{S} \neq λ_{S}, V_{0} K \cdot (\frac{1}{m t β_{S} V_{0}^{m} + K^{m}})^{\frac{1}{m}} & if β_{S} = λ_{S} . \end{matrix}

We note that in the extreme situation of having no nutrient supply (i.e. $S_{0} = 0$ ), it holds $δ^{+} (S_{0}) = 0$ and $δ^{-} (S_{0}) \approx 1$ . This results in $β_{S} = 0$ and $λ_{S} \approx λ + λ_{st} > 0$ . By inserting this special case of $β_{S} \neq λ_{S}$ into the given analytical solution, it simplifies to

V (t) S_{0} \to 0 ⟶ V_{0} K \cdot {⎛ ⎝ \frac{- (λ + λ_{st})}{- (λ + λ_{st}) K^{m} e^{(λ + λ_{st}) m t}} ⎞ ⎠}^{\frac{1}{m}} = V_{0} e^{- (λ + λ_{st}) t} .

This is an expected result, being the solution of a simple exponential decay model without growth:

{\begin{matrix} ˙ V & = - (λ + λ_{st}) V, V (0) & = V_{0} . \end{matrix}

Advantages of using the environmental stress level (ESL) instead to directly affect the cells.

For an alternative way to include the effect of the nutrient conditions in model $M_{opt}$ , we introduce an auxiliary variable: the environmental stress level $η = η (t)$ . It is designed to be an immeasurable quantity, describing to what extent viable cells would get stressed out, if exposed to certain environmental conditions, like in this case the nutrient saturation. We assume that the stress level is limited by a maximal level, which is set to one: $η \in [0, 1]$ .

In a general setting, the the ESL can be influenced by multiple environmental factors (e.g nutrient/oxygen/drug concentration), where each can be mathematically represented by a system variable. Let $E_{1} (t), \dots, E_{n} (t)$ be a set of such environmental variables, with

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ \begin{matrix} {˙ E}_{1} & = g_{1} (E_{1}, \dots, E_{n}, V, t), E_{1} (0) = E_{1, 0}, ⋮ {˙ E}_{n} & = g_{n} (E_{1}, \dots, E_{n}, V, t), E_{n} (0) = E_{n, 0}, \end{matrix}

being the given reaction equations. Then the ODE for the environmental stress level can be given by

˙ η = (n \sum j = 1 α_{j}^{-} \cdot δ_{j}^{-} (E_{j})) (1 - η)      \begin{matrix} increasing stress level (stressful conditions) \end{matrix} - (n \sum j = 1 α_{j}^{+} \cdot δ_{j}^{+} (E_{j})) η      \begin{matrix} recovery from stress (beneficial conditions) \end{matrix} with η (0) = η_{0} \in [0, 1] .

(2.3)

The first term includes the dynamics raising the ESL. Since $η \in [0, 1]$ , the increase is bounded by one from above via the factor $(1 - η)$ . If the environmental factors generate beneficial growth conditions for the cells, the ESL decreases according to the second term. The parameters $α_{j}^{-}$ resp. $α_{j}^{+}$ are variable-specific “sensitivity rates”, denoting how fast the stress level increases/decreases, if the value of the associated variable $E_{j}$ enters/leaves a range, which is critical for survival. With distinguishing between $α^{-}$ and $α^{+}$ per variable, it is possible to mathematically consider that cells react e.g. faster to the lack of an environmental factor, which is essential for their survival, than they recover from the deprivation once beneficial conditions are restored. Which values of $E_{j}$ are considered to be critical, is given by the explicit shape of the corresponding “influence functions” $δ_{j}^{-} (E_{j})$ and $δ_{j}^{+} (E_{j})$ . In particular, large values of $δ_{j}^{-} (E_{j})$ promote stress, whereas large values of $δ_{j}^{+} (E_{j})$ inhibits stress.

Collecting the influences of the environment in the stress level gives rise to some practical advantages. We can choose an arbitrary number of environmental variables and include them easily into the stress reaction due to the modular structure of (2.3). With the parameters $α_{j}^{-}$ and $α_{j}^{+}$ we can consider different sensitivities of the cells to changes in different environmental factors. Although neither these parameters nor the stress level $η$ can be measured directly in experiments, appropriate data of cell viability can provide enough information to estimate the values of the sensitivity rates. For this e.g. statistical methods for parameter calibration can be used, where the only prior assumption is a rough range, in which the parameter values are expected. Then, by grouping these sensitivity rates according to their magnitude, we can distinguish the corresponding environmental variables between having a fast or slow or even no impact on the cells. Under a quasi-steady state assumption, the ODE for $η$ can then be simplified by choosing a time scale of interest and omitting environmental variables with no influence.

In the considered experimental setting, we have the special case of the constant nutrient saturation $S (t)$ being the only environmental variable:

E_{1} (t) = S (t) = S_{0} with g_{1} (E_{1}, V, t) = g (S, V, t) = 0 .

We assume the ESL to increase and decrease with the same “nutrient sensitivity rate” $α_{S}$ . The corresponding “nutrient influence functions” are then $δ^{-} (S)$ and $δ^{+} (S)$ , as introduced previously in (2.2) for model $M_{S}$ . Hence, we have

α_{1}^{-} = α_{1}^{+} = α_{S}, δ_{1}^{-} (E_{1}) = δ^{-} (S_{0}) and δ_{1}^{+} (E_{1}) = δ^{+} (S_{0}) = 1 - δ^{-} (S_{0}),

i.e. given the initial condition $η (0) = η_{0} \in [0, 1]$ we consider the ODE

˙ η = α_{S} δ^{-} (S_{0}) \cdot (1 - η) - α_{S} (1 - δ^{-} (S_{0})) \cdot η = α_{S} (δ^{-} (S_{0}) - η) .

Since this equation is independent of the variable $V$ , the exact solution can be calculated by separation of variables and for $t \geq 0$ it is given as:

η (t) = δ^{-} (S_{0}) \cdot (1 - e^{- α_{S} t}) + η_{0} e^{- α_{S} t} .

(2.4)

Instead of the nutrient-dependent scaling functions $δ^{\pm}$ in the basic model $M_{S}$ , now the ESL influences the proliferation and starvation rate. A high stress level causes slower proliferation and faster starvation, whereas a low stress level has the opposite effect. Since by definition $η (t) \in [0, 1]$ , this results in the following initial value problem:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} ˙ V & = (1 - η (t)) \cdot β      = β_{η} (t) V (1 - {(\frac{V}{K})}^{m}) - (λ + η (t) \cdot λ_{st})      = λ_{η} (t) V, ˙ η & = α_{S} δ^{-} (S_{0}) \cdot (1 - η) - α_{S} (1 - δ^{-} (S_{0})) \cdot η, V (0) & = V_{0}, η (0) = η_{0} . \end{matrix}

(

M_{η}

)

By inserting the explicit time-dependent analytical solution (2.4), the system reduces to a non-autonomous ODE for $V$ , which is solved numerically. We observe that model $M_{η}$ has an interesting relation to the model $M_{S}$ : under the assumption that critical nutrient changes immediately influence the viable cells, i.e. assuming $α_{S} \to \infty$ , we return to model $M_{S}$ , since for $t \geq 0$ it holds

η (t) (???) = δ^{-} (S_{0}) \cdot (1 - e^{- α_{S} t}) + η_{0} e^{- α_{S} t} α_{S} \to \infty ⟶ δ^{-} (S_{0}) = const. \Rightarrow β_{η} (t) α_{S} \to \infty ⟶ β_{S}, λ_{η} (t) α_{S} \to \infty ⟶ λ_{S} .

The same result can be achieved by a quasi-steady state assumption, stating that changes in the ESL happen at a much faster time scale than changes in the tumor cell number. In this case, the stress level reaches its steady state virtually instantly, which is:

˙ η = α_{S} (δ^{-} (S_{0}) - η) = 0 \Rightarrow η = δ^{-} (S_{0}) .

We are aware that having multiple environmental factors included into the model would make the advantages of using the concept of the environmental stress level more evident. However, due to the lack of corresponding data, we focus on the special case of having only one environmental variable (nutrient saturation $S$ ). We consider the investigations in this manuscript as a proof of principle that the environmental stress level is a feasible alternative way to model the effect of the environment on the cells.

Overview over all models and their mathematical properties.

The last model $M_{η}$ has been constructed from model $M_{S}$ , which itself uses model $M_{opt}$ as a basis. Therefore, the number of variables and parameters increases with adding more complex dynamics. An overview over the variables and strictly positive model parameters of each system as well as abbreviating notations for important parameter terms/functions can be found in Tables 1, 2 and 3 below. The units of the variables and parameters are motivated by the biological setting and the measurement methods, which are explained in more detail in the following Section 2.2.

$M_{opt}$	$M_{S}$	$M_{η}$	Variable	Meaning	Unit
✓	✓	✓	$V = V (t)$	Density of viable tumor cells	$$ 10^{5} $ cells / m l$
	✓	✓	$S = S (t)$	Nutrient saturation	$\frac{$ 10 $ % FBS}{}$
		✓	$η = η (t)$	Environmental stress level	–

Table 1: Overview over the model variables and their occurrence in the models.

$M_{opt}$	$M_{S}$	$M_{η}$	Parameter	Meaning	Unit
✓	✓	✓	$β$	Maximal possible proliferation rate	$1 / day$
✓	✓	✓	$K$	Carrying capacity of the biological system	$$ 10^{5} $ cells / m l$
✓	✓	✓	$1 / m$	Strength of proliferation contact inhibition	–
✓	✓	✓	$λ$	Natural death rate	$1 / day$
	✓	✓	$λ_{st}$	Maximal possible starvation rate	$1 / day$
	✓	✓	$S_{thr}$	Nutrient sensitivity threshold	$\frac{$ 10 $ % FBS}{}$
		✓	$α_{S}$	Sensitivity rate of nutrient changes on stress level	$1 / day$

Table 2: Model parameters, their occurrence in the models and their meaning.

$M_{S}$	$M_{η}$	Notation		Meaning
✓	✓	$δ^{+} (S_{0}) =$	$\frac{S_{0}^{2}}{S_{thr}^{2} + S_{0}^{2}}$	Influence function for nutrient-promoted dynamics
✓	✓	$δ^{-} (S_{0}) =$	$1 - δ^{+} (S_{0})$	Influence function for nutrient-inhibited dynamics
✓		$β_{S} =$	$δ^{+} (S_{0}) \cdot β$	Nutrient-dependent net proliferation rate (constant)
✓		$λ_{S} =$	$λ + δ^{-} (S_{0}) \cdot λ_{st}$	Nutrient-dependent net death rate (constant)
	✓	$β_{η} (t) =$	$(1 - η (t)) \cdot β$	Stress-dependent net proliferation rate (varying in time)
	✓	$λ_{η} (t) =$	$λ + η (t) \cdot λ_{st}$	Stress-dependent net death rate (varying in time)

Table 3: Abbreviating notations for important/useful parameter terms and functions as well as their occurrence in the model equations (none of them occurs in model

M_{opt}

A mathematical analysis of the models yields positivity and boundedness of the solutions, which are important features for biological reasonableness. Table 4 summarizes the computed bounds as well as the steady states $¯ V$ reps. $(¯ V, ¯ η)^{T}$ and their stability. The corresponding calculations providing these results can be found in the supplementary Section 5.2.

Model		of the solutions	stable	unstable
		Positivity & boundedness		Steady states
$M_{opt}$		$0 \leq V (t) \leq max {V_{0}, K \sqrt[m]{1 - \frac{λ_{S}}{β_{S}}}}$	$K \sqrt[m]{1 - \frac{λ_{S}}{β_{S}}}$	$0$
$M_{S}$	for $λ_{S} < β_{S}$ :	$0 \leq V (t) \leq max {V_{0}, K \sqrt[m]{1 - \frac{λ_{S}}{β_{S}}}}$	$K \sqrt[m]{1 - \frac{λ_{S}}{β_{S}}}$	$0$
	for $λ_{S} \geq β_{S}$ :	$0 \leq V (t) \leq V_{0}$	$0$
$M_{η}$	for all $λ_{S}, β_{S}$ :	$0 \leq η (t) \leq max {η_{0}, δ^{-} (S_{0})} \leq 1$
	for $λ_{S} < β_{S}$ :	$0 \leq V (t) \leq max {V_{0}, K \sqrt[m]{1 - \frac{λ_{S}}{β_{S}}}}$	$⎛ ⎜ ⎝ \begin{matrix} K \sqrt[m]{1 - \frac{λ_{S}}{β_{S}}} δ^{-} (S_{0}) \end{matrix} ⎞ ⎟ ⎠$	$(\begin{matrix} 0 δ^{-} (S_{0}) \end{matrix})$
	for $λ_{S} \geq β_{S}$ :	$0 \leq V (t) \leq V_{0}$	$(\begin{matrix} 0 δ^{-} (S_{0}) \end{matrix})$

Table 4: Bounds of the solutions

V (t)

η (t)

and steady states

¯ V

(¯ V, ¯ η)^{T}

of each ODE model.

2.2 Experimental data

We use data from the experiments described in [7] to calibrate the unknown parameters from Table 2. A CellTiter-Blue^® assay was used to monitor the viability of tumor cells. In particular, viable cells metabolize a provided chemical and they emit measurable light as a result of this process. Hence, the data points are fluorescence intensity measurements. This way, viability was measured once every day and there are four biological replicates of each measurement to estimate statistical significance and repeatability.

Let $I^{total}$ be an intensity measurement of a specific cell line at time point $~ t$ . Excluding the corresponding background intensity $I^{BG}$ of the cell-free medium, the fluorescence intensity produced by viable cells $I^{V}$ is assumed to be directly proportional to the density of viable tumor cells $V (~ t)$ . The experiments were performed with different initial cell densities between $10^{4}$ and $10^{5}$ cells per milliliter, which motivates the unit of $V$ and eventually leads to the relation

I^{V} = I^{total} - I^{BG} \propto V (~ t) \frac{$ 10^{5} $ cells}{m l} \Rightarrow n = \frac{I^{V}}{V (~ t)} \in (0, 1),

where $n$ denotes the proportionality constant translating fluorescence intensity to cell density. Whenever we refer to “intensity measurements” in the following, we mean the fluorescence produced by the cells $I^{V}$ and neglect the superscript “ $V$ ” for better readability.

For nutrition the cells are supplemented with a particular concentration between $0 %$ to $10 %$ of fetal bovine serum (FBS). A supplementation with $0 \frac{% FBS}{}$ does not provide the cells with any nutrients, whereas $10 \frac{% FBS}{}$ generates optimal growth conditions. The nutrient supply is kept constant throughout the whole duration of each experiment, i.e. the nutrient saturation is assumed to stay at its initial level $S_{0} \in [0, 1]$ at any time. Overall, this motivates:

S (t) = S (0) = S_{0} \forall t \geq 0 with unit [S] = \frac{$ 10 $ % FBS}{} .

The remaining variable, the environmental stress level $η = η (t)$ , is an auxiliary variable and especially immeasurable. Hence, it has no experimental counterpart and is used as a dimensionless quantity.

How are uncertainties considered?

In reality, the equation $n = I / V (~ t)$ is not rigorously fulfilled. This can be due to e.g. model inadequacy or biological fluctuations of the cells’ metabolism, which affect measurement accuracy. To capture this uncertainty, we assume a multiplicative noise for each element of a set of $M \in N$ measurements ${I_{i}}_{i = 1}^{M}$ and a set of model solutions ${V_{i}}_{i = 1}^{M}$ considering the corresponding values of $t$ , $V_{0}$ , and $S_{0}$ . In the Bayesian framework that we adopt (see Section 2.3

) we model this to be a random variable

ε_{i}

such that for

i = 1, \dots, M

I_{i} = n V_{i} \cdot ε_{i} .

(2.5)

Let $ε_{i}$ be i.i.d. and have the unimodal and continuous distribution of a random variable $ε$

with probability density function (PDF)

f_{ε}

. Then, the following properties should hold:

	$[P1] :$	$supp (f_{ε}) \subseteq R^{+}, i.e. % all measurements are positive;$
	$:$	$E [ε] = 1, i.e. measurements are accurate % on average;$
	$:$	$f_{ε} (x) x \to \infty ⟶ 0 % and f_{ε} (x) x \to 0 ⟶ 0, i.e. % outliers are possible but not likely .$

Different distributions are possible to accomplish these properties. A small number of shape parameters and an easy calculation to ensure property [P2] are desirable. Therefore, we choose a Gamma distribution

ε \sim Γ (a, b)

a, b > 0

with a few restrictions. This is a plausible choice for multiplicative noise and often used in imaging theory (see e.g. [21, 22, 23]). Property [P1] and the desired behavior near infinity of [P3] are satisfied by definition. The corresponding PDF is given by

f_{ε}^{a, b} (x) = \frac{b^{a}}{Γ (a)} x^{a - 1} e^{- b x}

, where

Γ (\cdot)

denotes the Gamma function. To fulfill the remaining properties, we observe

E[ε]=ab\tiny[P2]=1, if %  b=a⇒ limx→0fa,aε(x)=limx→0(aaΓ(a)xa−1e−ax)\tiny[P3]=0, if  a>1.

Note that by constraining $b = a > 1$ the shape of the distribution depends only on the parameter $a$ . In fact, $a$

is directly related to the standard deviation

σ

and hence the variance of the distribution:

σ^{2} = Var (ε) = \frac{a}{a^{2}} = \frac{1}{a} \Rightarrow a = \frac{1}{σ^{2}}

. Therefore, for

i = 1, \dots, M

the uncertainty factor

ε_{i}

for a particular measurement

I_{i}

can be modeled by

ε_{i} = \frac{I_{i}}{n \cdot V_{i}} \sim Γ (\frac{1}{σ^{2}}, \frac{1}{σ^{2}}) with σ^{2} \in (0, 1) .

(2.6)

We use the percentiles $P_{5 %}$ and $P_{95 %}$ of the Gamma distributed uncertainty factors $ε_{i}$ to define the

``90\% (uncertainty) range'' around the solution V:[V⋅P5%,V⋅P95%].

(2.7)

In particular, given a specific noise variance $σ^{2}$ , the model expects 90% of the measurements within this interval, whereas respectively 5% are expected below and above it. The left side of Figure 2 depicts exemplary plots of the PDF $f_{ε}$ for different values of $σ^{2}$

. It also shows the positive skewness of the distribution. Under consideration of the measurement method, this is a reasonable feature for the uncertainty factors, assuming the cells might not metabolize the assay to their full potential. On the right side of Figure

2 we see an example of a 90% uncertainty range around a solution.

Figure 2: Left: Probability density function $f_{ε}$ of $ε \sim Γ (1 / σ^{2}, 1 / σ^{2})$ for varying variances $σ^{2}$ and the percentiles $P_{5 %}$ and $P_{95 %}$ for $σ^{2} = 0.05$ . Right: Resulting 90% uncertainty range (shaded area) around an exemplary model solution $V (t)$

The reason for considering a multiplicative (and, for instance, not additive) noise term is twofold. From a mathematical perspective, it allows to preserve positivity of the data. A more practical motivation is the reasonable assumption that the fluorescence noise of the intensity measurements is proportional to the density of viable cells. This is also supported by the observation of a larger variance for experiments with larger cell numbers in our data. Furthermore, multiplicative noise has demonstrated to be better suited for fluorescence than additive noise in other experimental settings [24].

How can the measurements be utilized?

One of the experiments in [7] monitors the dependence between cell viability and nutrient supply by providing five different nutrient concentrations $S_{0}$ over a period of 7 days: this results in five data sets “D1”–“D5”. These measurements are employed to calibrate the unknown model parameters for models $M_{S}$ and $M_{η}$ with Bayesian inversion methods (for details see following Section 2.3). Another experiment in [7] was designed to investigate cell behavior under optimal growth conditions, i.e. the cells were provided with $10 \frac{% FBS}{}$ over a period of 21 days. We denote the data set generated by this experiment with “D6”. It is used additionally to D1–D5 to validate the calibration results with model $M_{opt}$ . All experiments start with several populations of different initial size $V_{0}$ . Table 5 shows the experimental setup of all data sets and the resulting initial conditions in the models.

Data set	(in $\frac{% FBS}{}$ )	(in days)	$V_{0}$	$S_{0}$	$η_{0}$
	Nutrition	Duration	Initial values in the models
D1	10.0	7	$1.00, 0.50, 0.25$	$1.00$	$0.00$
D2	07.5	7	$1.00, 0.50, 0.25$	$0.75$	$0.00$
D3	05.0	7	$1.00, 0.50, 0.25$	$0.50$	$0.00$
D4	02.5	7	$1.00, 0.50, 0.25$	$0.25$	$0.00$
D5	00.0	7	$1.00, 0.50, 0.25$	$0.00$	$0.00$
D6	10.0	21	$1.00, 0.50, 0.25, 0.10, 0.05$	( $1.00$ )	( $0.00$ )

Table 5: Overview over the data sets and the corresponding initial values. D1–D5 are used to calibrate the models

M_{S}

and

M_{η}

and D6 for validation of the results with model

M_{opt}

2.3 Parameter calibration: Bayesian inversion and Sequential Monte Carlo (SMC)

The task to identify the unknown true parameters from given data is called the “inverse problem”. We solve it using a Bayesian approach, which leads, under mild assumptions, to a naturally well-posed inverse problem [25]. Furthermore, this approach allows to quantify the uncertainty in the estimated parameters and hence it brings more information than a deterministic inversion method, at the price of a higher computational cost.

We collect all parameters to be estimated in a vector

θ \in Θ \subseteq R^{d}

(

d \in N

), where

Θ

denotes the parameter space. Assuming the considered data consist of

M

intensity measurements (excluding background intensity), we collect them in a vector . Defining

G : Θ \to R^{M}

with

θ \mapsto (G_{i} (θ))_{i = 1}^{M}

as the forward operator mapping parameter values to the corresponding intensities, such a measurement can be rewritten as

I_{i} (???) = n V_{i} \cdot ε_{i} = G_{i} (θ) \cdot ε_{i},

(2.8)

where $V_{i}$ is the corresponding model solution to ( $M_{opt}$ ), ( $M_{S}$ ) or ( $M_{η}$ ) using $θ$ as parameters and $ε_{i}$ is the multiplicative noise. Note that we use parts of the data sets D1–D5 for parameter estimation, hence each measurement $I_{i}$ can refer to a different nutrient condition, initial cell density, and time point. Therefore, $V_{i}$ has to be calculated in consideration of the corresponding values for $S_{0}$ , $V_{0}$ and $t$ .

Now, we want to consider all measurements in $I = (I_{i})_{i = 1}^{M}$ collectively. The parameter vector $θ$ and the noise $(ε_{i})_{i = 1}^{M}$ are modeled as multi-dimensional random variables taking values in $Θ$ and $R^{M}$ respectively. The Bayesian formulation of the problem is the following: Given a prior (measure) $μ_{0}$ on $Θ$ , compute the posterior (measure) $μ^{I}$ given the data $I$ . The prior in the Bayesian setting is the correspondent of a regularization in deterministic inverse problems [25] and it reflects the knowledge about the parameters before including any information given by the data, whereas the posterior describes the knowledge after seeing the data. Let $π_{0}$ and $π^{I}$ denote the probability densities of $μ_{0}$ and $μ^{I}$ , respectively. By Bayes’ formula, we have

π^{I} (θ) = \frac{L (I | θ) \cdot π_{0} (θ)}{} \int_{Θ} L (I | θ) \cdot π_{0} (θ) d θ \propto L (I | θ) π_{0} (θ),

(2.9)

where $L$ is the data likelihood [25]. The proportionality constant of relation (2.9) depends only on $I$ . It is called “model evidence” and it can be used to quantitatively compare two models (see Section 2.4). We remind from the previous Section 2.2 that (2.6) states i.i.d. $ε_{i} \sim Γ (\frac{1}{σ^{2}}, \frac{1}{σ^{2}})$ with $σ^{2} \in (0, 1)$ for every measurement $I_{i}$ , $i = 1, \dots, M$ . Using this together with (2.8), the data likelihood of $I = (I_{i})_{i = 1}^{M}$ is:

L (I | θ) \propto M \prod i = 1 {(\frac{I_{i}}{G_{i} (θ)})}^{1 / σ^{2} - 1} exp (- \frac{1}{σ^{2}} \frac{I_{i}}{G_{i} (θ)}) .

How to sample from the posterior?

In order to make predictions, we want to sample from the given posterior distribution (2.9). However, this has a complicated, concentrated density, so we cannot sample from it exactly with a random number generator. To approximate the posterior measure, we use therefore the Sequential Monte Carlo (SMC) method, which we explain now based on [26]. In SMC one considers a sequence of intermediate distributions $(μ_{k})_{k = 0}^{N}$ , such that $μ_{0}$ is the prior and $μ_{N} = μ^{I}$ coincides with the posterior distribution. The probability density $π_{k}$ of the intermediate measure $μ_{k}$ can be defined by

π_{k} (θ) = \frac{1}{Z_{k}} k \prod i = 1 L_{i} (I_{i} | θ) π_{0} (θ) or, equivalently, π_{k} (θ) = \frac{1}{{~ Z}_{k}} L_{k} (I_{k} | θ) π_{k - 1} (θ),

(2.10)

where $Z_{k}$ and ${~ Z}_{k}$ are normalizing constants and

L_{i} (I_{i} | θ) = \frac{1}{Z_{i}} {(\frac{I_{i}}{G_{i} (θ)})}^{1 / σ^{2} - 1} exp (- \frac{1}{σ^{2}} \frac{I_{i}}{G_{i} (θ)})

is the likelihood associated to observation $I_{i}$ with normalization constant $Z_{i}$ , $i = 1, \dots, M$ . The intermediate densities could also be constructed with an adaptive approach using tempering [9]. However, since the considered data measures the quantity of interest in a time series, the presented filtering method is computationally more efficient. The SMC algorithm samples sequentially from the intermediate measures $μ_{k}$ using a weighted swarm of samples, called particles. Let $P \in N$ be the sample size, i.e. the number of particles. At the $k$ -th iteration ( $k = 1, \dots, N$ ) the algorithm leads to a collection of particles ${θ_{p}}_{p = 1}^{P}$ with associated weights ${W_{p}^{k}}_{p = 1}^{P}$ , which gives the approximation

π_{k} (θ) \approx P \sum p = 1 W_{p}^{k} δ_{θ_{p}} (θ) .

The SMC algorithm is summarized in Algorithm 1: we achieve appropriately weighted particles to approximate $μ_{N} = μ^{I}$ by starting with uniformly weighted particles distributed according to the prior $μ_{0}$ (line 1) and iteratively move the samples from the previous measure $μ_{k - 1}$ to $μ_{k}$ in a selection (lines 3-7) and a mutation step (line 8) [27, Ch. 5], which are explained in detail in the following paragraphs.

k = 0

: sample

θ_{p} \sim μ_{0}

and set

W_{p}^{0} = \frac{1}{P}, p = 1, \dots, P

2:for

k = 1, \dots, N

w_{p}^{k} = L_{k} (I_{k} | θ) W_{p}^{k - 1}, W_{p}^{k} = w_{p}^{k} / (\sum_{p = 1}^{P} w_{p}^{k}), p = 1, \dots, P

▹

selection step

4: if

P_{eff} < P_{thr}

then resample:

(i)

: sample indices acc. to distribution

R

of particle indices:

(π_{1}, \dots, π_{P}) \sim R (W_{1}^{k}, \dots, W_{P}^{k})

(i i)

: set

θ_{p} = θ_{π_{p}}

and

W_{p}^{k} = \frac{1}{P}, p = 1, \dots, P

7: end if

8: move

θ_{p} \sim κ_{k} (θ_{p}, \cdot), p = 1, \dots, P

▹

μ_{k}

-invariant mutation step

9:end for

Algorithm 1 Sequential Monte Carlo

Selection step. We start with a collection of particles ${θ_{p}}_{p = 1}^{P}$ , distributed according to $μ_{k - 1}$ . Their weights ${W_{p}^{k - 1}}_{p = 1}^{P}$ are updated to ${W_{p}^{k}}_{p = 1}^{P}$ by importance sampling: for $p = 1, \dots, P$ we have

W_{p}^{k} = \frac{w_{p}^{k}}{\sum_{p = 1}^{P} w_{p}^{k}} with w_{p}^{k} (???) = L_{k} (I_{k} | θ) W_{p}^{k - 1} .

We see that $W_{p}^{k}$ is normalized to ensure that $\sum_{p = 1}^{P} W_{p}^{k} = 1$

, i.e. having a probability distribution. Note that the importance sampling only changes the weights and not the particles. However, if there are many particles with low weights, the estimation is only as accurate as a Monte Carlo approximation with a very small number of particles

[28]. In this case, the reweighing step is followed by a resampling step, where the particles are replaced according to their updated weights. Resampling is needed, if the effective sample size

P_{eff} \approx {(P \sum p = 1 {(W_{p}^{k})}^{2})}^{- 1}

is small, which we check by comparison with a threshold $P_{thr} = τ P, τ \in (0, 1)$ :

This discards particles with low weight and improves the representation of the distribution $μ_{k}$ .

Mutation step. Performing only selection steps will eventually lead to degeneracy in the diversity of the particle population. In particular, after some resampling steps, few particles will survive and be replicated. Therefore, we introduce diversity in the particles by moving them according to a Markov Chain Monte Carlo (MCMC) kernel $κ_{k} (\cdot, \cdot)$ . This kernel is $μ_{k}$ -invariant, i.e. it does not modify the particle distribution. We adopt the adaptive strategy developed in [9] to construct such a MCMC kernel. A random walk Metropolis-Hastings (MH) proposal is used on each univariate component, conditionally independently. More precisely, remembering that each particle is a vector $θ_{p} \in R^{d}$ , MH proposes ${q_{p}}_{p = 1}^{P}$ with $q_{p} \in R^{d}$ by computing

(q_{p})_{j} = (θ_{p})_{j} + ϵ_{j}^{k} \cdot ξ_{j} with ξ_{j} \sim N (0, 1),

where $(\cdot)_{j}$ denotes the $j$ -th component ( $j = 1, \dots, d$ ) and the scale $ϵ_{j}^{k}$ is tuned to the acceptance rate of the previous SMC iteration. This means, we choose $ϵ_{j}^{k}$ adaptively as

ϵ

Environmental stress level to model tumor cell growth and survival

Modeling glioma invasion with anisotropy- and hypoxia-triggered motility enhancement: from subcellular dynamics to macroscopic PDEs with multiple taxis

Evaluation of Entropy and Fractal Dimension as Biomarkers for Tumor Growth and Treatment Response using Cellular Automata

Holographic Cell Stiffness Mapping Using Acoustic Stimulation

Parallel Model Exploration for Tumor Treatment Simulations

What do cells regulate in soft tissues on short time scales?

Modeling Cell Populations Measured By Flow Cytometry With Covariates Using Sparse Mixture of Regressions

Antioxidant capacity is repeatable across years but does not consistently correlate with a marker of peroxidation in a free-living passerine bird

1 Introduction

2 Materials and method

2.1 Mathematical modeling

What if the cells find optimal nutrient conditions?

What if suboptimal nutrient conditions directly affect the cell dynamics?

Advantages of using the environmental stress level (ESL) instead to directly affect the cells.

Overview over all models and their mathematical properties.

2.2 Experimental data

How are uncertainties considered?

How can the measurements be utilized?

2.3 Parameter calibration: Bayesian inversion and Sequential Monte Carlo (SMC)

How to sample from the posterior?

Environmental stress level to model tumor cell growth and survival

Related Research

Modeling glioma invasion with anisotropy- and hypoxia-triggered motility enhancement: from subcellular dynamics to macroscopic PDEs with multiple taxis

Evaluation of Entropy and Fractal Dimension as Biomarkers for Tumor Growth and Treatment Response using Cellular Automata

Holographic Cell Stiffness Mapping Using Acoustic Stimulation

Parallel Model Exploration for Tumor Treatment Simulations

What do cells regulate in soft tissues on short time scales?

Modeling Cell Populations Measured By Flow Cytometry With Covariates Using Sparse Mixture of Regressions

Antioxidant capacity is repeatable across years but does not consistently correlate with a marker of peroxidation in a free-living passerine bird

1 Introduction

2 Materials and method

2.1 Mathematical modeling

What if the cells find optimal nutrient conditions?

What if suboptimal nutrient conditions directly affect the cell dynamics?

Advantages of using the environmental stress level (ESL) instead to directly affect the cells.

Overview over all models and their mathematical properties.

2.2 Experimental data

How are uncertainties considered?

How can the measurements be utilized?

2.3 Parameter calibration: Bayesian inversion and Sequential Monte Carlo (SMC)

How to sample from the posterior?