1. Introduction

Industrial chemicals might be environmental or health hazards, thus there is an increasing need for developing efficient methodologies for toxicity assessment. The traditional method of conducting in vivo assays is very efficient for understanding the mode of action of various toxicants, but it is expensive, time consuming, and there are moral issues regarding the use of animals. Cell-based in vitro assays [28] are attractive alternatives for the assessment of toxicity. We propose a mathematical model to study the effect of toxicants on human cells.

Our work is part of the cytotoxicity profiling project carried by Alberta Centre for Toxicology in which initially 63 chemicals [27] were investigated using the xCELLigence Real-Time Cell Analysis High Troughput (RTCA HT) Assay [19]. The xCELLigence RTCA HT system is developed by the ACEA Biosciences Inc. (San Diego, USA), utilizing 384 well electronic microtiter plates (E-Plates 384). Impedance is generated by cells attached to electrodes, impeding the electric current between electrodes. Therefore, the impedance signal is label free and non-invasive, enabling direct measurement of cellular status in real time. The impedance value is converted by a software to Cell Index, which closely reflects not only cell number, but also cell morphology and adhesion.

The cells were seeded into the wells of E-Plates 384 [19], and the test substances were added about 24 hours later to allow cells to attach to the bottom of the wells and adapt to growth. The 11 concentrations of the test substances were selected based on a literature study, and each of them was tested in at least quadruplicates [27]. Time-dependent concentration response curves (TCRCs) for each test substance in each cell line were generated [27].

Series of physiological events, such as changes in cell population, cell morphology, and cellular functions, that characterize an adverse biological response are referred to as mode of action (MOA) [27]. According to the MOA, the tested substances were divided in 10 groups [27]. For example, the TCRCs for two chemicals from the same group, PF431396 and monastrol, are displayed in Fig. 1.

The main objectives of our work are to develop a mathematical model capable of reproducing the TCRCs and, based on this model, to investigate the effects of toxicants on the survival/death of the cells' population. We consider a three dimensional model for a cells' population subjected to an acute dose of toxicant. The state variables are the cell index, as a proxy for the number of cells, the concentration of toxicant in the cell, and the concentration of toxicant in the environment, and they are coupled by a linear dose-response function. This model is a special case of the class of models proposed in [4], and it is related to the models considered in [5,10,15]. However, since we consider an acute dose of toxicant instead of a chronic one, the analysis of the survival/death of the cells' population is different from the one done in the previously mentioned papers. We analyze the stability of this system of ordinary differential equations analytically using an approach similar with the one in [7], and we find sufficient conditions for the survival or death of the cells' population.

An important contribution of this paper is the efficient statistical methodology used to validate the proposed mathematical model. For parameter estimation and prediction we propose an approach based on the Expectation Maximization (EM) algorithm [3] and the Unscented Filter [12]. This method can also be used to predict the concentration of toxicant outside the cells. These predictions are important in any system for toxicity monitoring. For example, a similar cytotoxicity mathematical dynamic model [8] is applied for early detection and quantification of the presence of toxicants in the water supply system. The estimation of the model parameters and prediction of the concentrations of the toxic agent presented in water supply are done using the moving horizon estimator (MHE) and extended Kalman filer (EKF). Although the model considered in [8] is less complicated than the model presented in this paper, for parameter estimation the MHE and the EKF require extensive tunning, so the EM algorithm is a much more efficient alternative.

Since the parameters of the proposed model are estimated not measured directly, there is uncertainty regarding their precise values. We use Monte Carlo simulations and parameter sensitivity analysis for uncertainty quantification. Sensitivity analysis [13] refers to a broad group of methods that provide valuable insight regarding the dependence of the model output on its parameters. We apply a global sensitivity analysis (GSA) method [23] to test the robustness of the results of our model in the presence of uncertainty, and to increase the understanding of the relationships between input and output variables in the model.

In addition to proposing an efficient statistical method for parameter estimation and prediction of the concentration of toxicant outside the cells, an important application of the mathematical model is to generate numerically the TCRCs. Based on these numerical simulations and a mathematical analysis of the model, we can determine the threshold value for the initial concentration of a chemical after which the death of the cells' population happens (i.e. the chemical considered becomes toxic for the cells).

From the information included in the TCRCs, various parameters, such as LC$_{50}$ [14], KC$_{50}$ [29,20], and AUC$_{50}$ [19], were used to describe the time and concentration dependent cellular activity, and they were used in various statistical methods for clustering [27]. The success of clustering and classification methods, such as neural networks and support vectors machines, depends on providing TCRCs corresponding to a range of values of the initial concentration of the toxicant that include both cases when the concentration is low enough to not affect the cells' population, and cases when the concentration is high enough to determine the eventual death of the cells' population ([20,27]). For example, in Fig. 1(b) none of the concentrations used for monastrol was high enough to kill cells. Using the proposed mathematical model we can determine an appropriate range for the initial concentration of the toxicant. Our methodology based on stability analysis has better mathematical justification and gives more accurate results than the approach proposed in [20] and based on LC$_{50}$ and linear regression.

The paper is organized as follows. We start by presenting the mathematical model in section 2 followed by the study of stability in section 3. Section 4 includes the methods used for parameter estimation and prediction. Next, in section 5, we validate the proposed model by comparing the experimental TCRCs with the predicted values for the cell index. The GSA and Monte Carlo Simulations are presented in section 6. The last section of the paper includes a summary and discussion of the main results and intended future work.

---

2. The mathematical model

The data set used in this study is from the cytotoxicity test on human hepatocellular carcinoma cell line HepG2 [19].

Without exposure to the test substances, a typical cell growth curve exhibits four phases: the lag, the log, the plateau and the decline phases. The lag phase is the phase when the cells do not divide and adapt to the culture conditions. In the log phase the cells actively proliferate and have an exponential increase in the cell density. In the plateau phase, the cellular proliferation reaches a steady-state level. Finally, in the decline phase, the cell death predominates and the number of viable cells declines due to the natural path of the cell cycle and depletion of the nutrient supplements.

We focus on the modeling of the log phase. After seeding the cells in E-Plates for 24 hours, eleven concentrations of each substance were applied. Negative controls were presented in each E-Plate 384, as the cells in these wells were treated with only assay buffer for substance dilution, but no testing substances. The cells were monitored totally for 89 hours. The TCRCs are truncated such that only the data from the log phase are kept.

To model the negative control curves, we start with the well-known logistic model describing cell growth in non-toxic environment [24]:

$\frac{dn(t)}{dt} = \beta n(t)(1-\frac{n(t)}{K}),\label{model0}$

(1)

where $n(t)$ is the cell index at time $t$, which closely reflects the number of cells, $\beta$ denotes the cell growth rate, and $K$ is the capacity volume. The initial stage of growth is approximately exponential; then, as the population is closed to capacity $K$, the growth slows.

We consider the following model for cells growth/death in a toxic environment [4]. Let $C_0(t)$ and $CE(t)$ be the the intracellular and the extracellular concentrations of toxicant at time $t$, respectively. We suppose the death rate of cell is linearly dependent on the intracellular concentration of toxicant with coefficient $\alpha>0$, so we have

$\frac{dn(t)}{dt} = \beta n(t)(1-\frac{n(t)}{K})- \alpha C_0(t) n(t). \label{model1}$

(2)

We suppose that the only source of toxicant is the environment. Linear kinetics is known to be relevant for certain toxicants, but there is evidence that higher order formulations are also important [4]. In [6] the mechanisms of in vitro cytotoxicity are analyzed, and in the proposed models the toxicant's uptake rate is considered to be a combination of a linear diffusive component and a saturable, carrier-mediated component (Michaelis-Menten like kinetics). The rapid on-line estimation of toxicant concentration in water supply proposed in [8] is also based on these models. However, after estimating the parameters, for some toxicants it appears that the Michaelis-Menten uptake is not significant [6], so we ignore carrier-mediated transportation and we use linear kinetics [4,24]:

$\frac{dC_0(t)}{dt} = \lambda_1^2 CE(t) -\eta_1^2C_0(t), \label{model2}$

(3)

where $\lambda_1^2$ represents the uptake rate of the toxicant from environment, and $\eta_1^2$ is the toxicant input rate to the environment.

As in [4] we assume that the cells absorb part of the toxicant, and we express the rate of change of the total amount of toxicant in the environment as the sum of the rates of input to and losses from the environment. Supposing that the inputs to the environment are proportional with $C_0(t)n(t)$, and the losses from the environment to the population are proportional with $CE(t)n(t)$, the concentration of toxicant outside the cells changes over time according the following equation:

$\frac{dCE(t)}{dt} = \lambda_2^2C_0(t)n(t) -\eta_2^2CE(t)n(t) \label{model31}$

(4)

where $\lambda_2^2$ is the toxicant uptake rate from cells, and $\eta_2^2$ represents the losses rate of toxicants absorbed by cells. Coupling the equations (2) - (4), we obtain the following toxicity model

$\begin{split} \frac{dn(t)}{dt} &= \beta n(t)(1-\frac{n(t)}{K})- \alpha C_0(t) n(t) \\ \frac{dC_0(t)}{dt} &= \lambda_1^2 CE(t) -\eta_1^2C_0(t)\\ \frac{dCE(t)}{dt} &= \lambda_2^2C_0(t)n(t) -\eta_2^2CE(t)n(t). \end{split}\label{model}$

(5)

The variables and parameters in this model are listed in Table 1.

Without loss of generality we suppose that the initial time is $t_0=0$, and we consider the initial cell index $0<n(0)\le K$. We assume that the toxicants do not exist in the cells before experiments, so the initial concentration of internal toxicant $C_0(0 )=0$, and the initial concentration of toxicants outside of the cells $CE(0 )>0$ equals to the concentration of toxicant used in the experiments.

As we have mentioned before, the model (5) is a particular case of the class of models introduce in [4]. More recently, similar models were studied in [5,10,15]. However, in these papers sufficient conditions for persistence or extinction of the population are expressed in terms of the toxicant input rate $u(t)$, so they do not apply to our model because we consider an acute dose of toxicant ($u=0$), instead of a chronic one. In the next section we study the stability of the model (5) applying both techniques for studying the persistence as in [4,5], and methods used to analyze the existence and stability of equilibria.

---

3. Stability analysis

We assume that all parameters of model (5) have positive values.

---

3.1. Existence and stability of equilibria

Equilibria are very important in the study of the asymptotic behavior of a dynamical system [7]. To determine the equilibria, we set the vector field in (5) equal to $\mathbf{0}$, and we solve the equilibria equations. To study the stability of the equilibria, we analyze the eigenvalues of the Jacobian matrix $J$ of (5):

$J=\begin{bmatrix} \beta-2\beta n/K-\alpha C_0 &-\alpha n &0\\ 0&-\eta^2_1 &\lambda^2_1\\ \lambda^2_2 C_0- \eta^2_2 CE &\lambda^2_2n&-\eta^2_2n \end{bmatrix}.$

We have the following equilibria.

• The extinction equilibria: $\begin{bmatrix} 0, CE^*\frac{\lambda_1^2}{\eta_1^2},CE^* \end{bmatrix}$ where $CE^*>0$ is the concentration of toxicant outside the cells. For $\begin{bmatrix} 0, CE^*\frac{\lambda_1^2}{\eta_1^2},CE^* \end{bmatrix}$ the eigenvalues are

$\xi_1 = \beta-\alpha CE^*\frac{\eta_1^2}{\lambda_1^2} , \xi_2 = -\eta_1^2, \xi_3 = 0$

To determine the stability and the qualitative behavior in a neighborhood of these non-hyperbolic equilibria, we study the behavior on the local center manifold [21] (see appendix A for details regarding the center manifold). The analytical center manifold is given by the line $\{(n, C_0,CE): n=0, \lambda_1^2 CE -\eta_1^2C_0=0\}$. The extinction equilibria are unstable [7] if the first eigenvalue is positive, i.e. $0<CE^*< \frac{\beta\eta_1^2}{\alpha\lambda_1^2}$.

• The maximum capacity equilibrium: $\begin{bmatrix} K, 0,0 \end{bmatrix}$. For $\begin{bmatrix} K, 0,0 \end{bmatrix}$ the eigenvalues are

$\begin{split} \xi_1 &= -\beta <0\\ \xi_{2,3} &= -\frac{1}{2}\left( \eta_1^2+\eta_2^2K \pm\sqrt{(\eta_1^2+\eta_2^2K)^2+4(\lambda_1^2\lambda_2^2-\eta_1^2\eta_2^2)K} \right) \end{split}$

Notice that if $\lambda_1^2\lambda_2^2-\eta_1^2\eta_2^2<0$, then all the eigenvalues are negative, so $\begin{bmatrix} K, 0,0 \end{bmatrix}$ is locally asymptotically stable. If $\lambda_1^2\lambda_2^2-\eta_1^2\eta_2^2>0$ then $\xi_1<0$, $\xi_2>0$, $\xi_3<0$, so $\begin{bmatrix} K, 0,0 \end{bmatrix}$ is unstable. Finally, if $\lambda_1^2\lambda_2^2-\eta_1^2\eta_2^2=0$ then $\xi_1<0$, $\xi_2=0$, $\xi_3<0$, so $\begin{bmatrix} K, 0,0 \end{bmatrix}$ is a non-hyperbolic fixed point.

• The interior equilibrium: If $\lambda_1^2\lambda_2^2-\eta_1^2\eta_2^2 = 0$, then we have the interior equilibria $\begin{bmatrix} K(1-\frac{\alpha \lambda_1^2 CE }{\beta \eta_1^2}), \frac{\lambda_1^2 CE}{\eta_1^2},CE \end{bmatrix}$, where $CE>\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$. For these equilibria the eigenvalues are

$\xi_1 = 0, \xi_2 = -\beta\left(1-\frac{\alpha \lambda_1^2 CE }{\beta \eta_1^2}\right)<0, \xi_3 = -\eta_1^2-\eta_2^2K\left(1-\frac{\alpha \lambda_1^2 CE }{\beta \eta_1^2}\right)<0,$

so these equilibria are non-hyperbolic fixed points.

---

3.2. Persistence and extinction

Following [5,4] we define persistence and extinction as follows.

• The cell population is said to be uniformly persistent if there exist the constants $M_1>0$, $M_2>0$ such that $0<M_1\le\lim\inf_{t\rightarrow\infty}n(t)\le \lim\sup_{t\rightarrow\infty}n(t)\le M_2<+\infty$.

• The cell population is said to go to local extinction if $\lim_{t\rightarrow \infty}n(t)=0$.

Intuitively, uniform persistence corresponds to the long term cells survival, and the cells death corresponds in the previous definition to local extinction.

Similarly with [5] we can prove the following lemma (see appendix B for a detailed proof).

Lemma 3.1. If $\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2>0$, then $0<CE(t)\le CE(0)$, $0\le C_0(t)\le \frac{\lambda_1^2CE(0)}{\eta_1^2}$, and $n(t)>0$, for all $t\ge 0.$

We have noticed that for the chemicals considered in this paper, the estimated values of the parameters $\eta_1$, $\eta_2$, $\lambda_1$, and $\lambda_2$ verify $\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2>0$, so we focus on this case. Using lemma 3.1, we get the following result (see appendix B for the proof).

Theorem 3.2. If $\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2>0$, then $\lim_{t\rightarrow \infty} CE(t)$ exists and its value determines the asymptotic behavior of system (5) according to the following two cases:

1. If $\lim_{t\rightarrow \infty} CE(t)<\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$ then the population is uniformly persistent and we have

$\lim\limits_{t\rightarrow\infty} n(t)=K, \lim\limits_{t\rightarrow\infty} C_0(t)=\lim\limits_{t\rightarrow\infty} CE(t)=0.$

2. If $\lim_{t\rightarrow \infty} CE(t)>\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$ then the population goes to local extinction and we have

$\lim\limits_{t\rightarrow\infty} n(t)=0, \lim\limits_{t\rightarrow\infty} C_0(t)=CE^*\frac{\lambda_1^2}{\eta_1^2}, \lim\limits_{t\rightarrow\infty} CE(t)=CE^*>\frac{\beta\eta_1^2}{\alpha\lambda_1^2},$

Moreover, the population component is in $L^1[0,\infty)$, i.e. $|n|_1=\int_0^\infty n(t)dt<\infty$.

Notice that if $\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2>0$ and $0<CE(0)<\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$ then, by Lemma 3.1, $\lim_{t\rightarrow\infty}CE(t)< \frac{\beta\eta_1^2}{\alpha\lambda_1^2}$, so we are in the first case and we have $\lim_{t\rightarrow\infty} n(t)=K>0$. From the analysis of the stability of the equilibria we know that $[K, 0, 0]$ is a locally asymptotically stable fixed point, so there exists populations modeled by equations (5) that are uniformly persistent. These results are summarized in the following corollary.

Corollary 1. If $\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2>0$ and $0<CE(0)<\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$ then the population is uniformly persistent, and we have $\lim_{t\rightarrow\infty}n(t)=K$.

This is illustrated numerically by the trajectories displayed in Fig. 2(a) for monastrol. These trajectories correspond to several sets of initial values $n(0)>0$, $C_0(0)=0$, $0<CE(0)<\frac{\beta\eta_1^2}{\alpha\lambda_1^2}=6.51$. We noticed that for this set of initial values the trajectories approach the equilibrium $[K,0,0]=[18.17,0,0]$.

This means that if the initial concentration is sufficiently low, the presence of the chemical does not affect the convergence of the cell index towards the maximum capacity $K$. Moreover, we have a threshold, $\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$, such that if the experiment is done using a chemical with initial concentration $CE(0)<\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$, then the cell population is uniformly persistent.

If $\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2>0$ and $CE(0)>\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$, then we know that the population goes to local extinction if $|n|_1<\infty$ (see Theorem 3.2). In Fig. 2(b) this is illustrated numerically for monastrol by the trajectories that asymptotically approach the extinction equilibria $\begin{bmatrix} 0, CE^*\frac{\lambda_1^2}{\eta_1^2},CE^* \end{bmatrix}$, with $CE^*> \frac{\beta\eta_1^2}{\alpha\lambda_1^2}$. We notice from Fig. 2(b) that these extinction equilibria are on a line representing also the center manifold.

However, if $\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2>0$, $CE(0)>\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$ and $|n|_1=\infty$, the population is uniformly persistent (see Theorem 3.2). Thus, the initial concentration of the chemical and the entire past history of the population are important to determine the persistence or extinction of the population. A similar conclusion was obtained in [4] for a reduced model (with $\lambda_2=0)$, also corresponding to an acute dose of toxicant.

For monastrol, the separation between persistence and extinction according to the initial values $n(0)>0$, $C_0(0)=0$, $CE(0)>\frac{\beta\eta_1^2}{\alpha\lambda_1^2}=6.51$ are presented in Fig. 3. Results of numerical simulations similar with the ones presented in Fig. 3 are important, because they can be applied to choose an appropriate range for the values of $CE(0)$ used in the experiments such that the TCRCs will illustrate the evolution from cells' population survival (corresponding to small values of $CE(0)$) to cells' population extinction (corresponding to large values of $CE(0)$).

In the next section we design an efficient method to estimate the parameters based on the TCRCs. With the proposed method, for a given chemical and a cell line, the parameters can be determined after only two experiments, namely the negative control and an experiment with a relatively high initial concentration of the testing substance. Thus, after two experiments we can determine the threshold $\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$ and we can obtain numerically separation regions similar with one displayed in Fig. 3. This information can be used to save both experimental time and resources.

---

4. Parameter estimation

The non-linear system of ordinary differential equations (5) cannot be solved analytically, so to estimate the parameters we rewrite it in a state-space form. Notice that from the experimental data recorded in the TCRCs we get observations, possibly affected by measurement errors, only for the cell index $n(t)$. Thus, using the Euler integration scheme with time step $h$, the associated discrete state-space system can be expressed as:

$x_{k+1}=x_k+h\begin{bmatrix} \beta x_k[1](1-\frac{x_k[1]}{K})- \alpha x_k[2] x_k[1] \\ \lambda_1^2 x_{k}[3] -\eta_1^2x_k[2] \\ \lambda_2^2x_k[2]x_k[1] -\eta_2^2x_k[3]x_k[1]\end{bmatrix} + v_{k+1}$

(6)

$y_{k+1}= C x_{k+1}+w_{k+1}$

(7)

Here $x_{k+1}=[n,C_0, CE]^{t}$ is the state of the system, $y_{k+1}$ is the observation at time step $k+1$, and $C=\begin{bmatrix} 1&0&0 \end{bmatrix}$. We assume that the white noise vectors, $v_k$ and $w_k$, are Gaussian and from uncorrelated white sequences such that $v_k \sim N(0,Q)$ and $w_k \sim N(0,R)$. The $3\times 3$ diagonal covariance matrix $Q$, the variance $R>0$, $\alpha$, $\lambda_1,\lambda_2, \eta_1$, $\eta_2$, $\beta$ and $K$ represent the unknown parameters, $h$ is the time step and the initial values $x_1[1]=n(0)$ comes from the measured data at the initial time $t_0$, $x_1[2]=C_0(0)=0$ and $x_1[3]=CE(0)$ equals the concentration of toxicant used in the experiments.

The parameters $\beta$ and $K$ are also the parameters of the logistic model (1) for which we can easily get an analytic solution. Thus, using the experimental data corresponding to the negative control (no toxicant) we can estimate $\beta$ and $K$ using the nonlinear least square method based on the analytic solution of equation (1).

The system (5) cannot be solved analytically, so to estimate the remaining parameters $\Theta=\{Q, R,\alpha, \lambda_1,\lambda_2, \eta_1, \eta_2\}$ we use the state-space model (6)-(7) and the Expectation Maximization (EM) algorithm based on the the unscented filter (UF). The EM algorithm is a classical method for estimating the parameters of linear systems [22], and for a general nonlinear dynamics, the EM algorithm has been applied using the UF or the EKF smoothing [3,25].

The EM algorithm is especially useful when it is easier to calculate the likelihood of the model using not only the observed data $Y_{obs}=\{y_k, k=1,\ldots,N\}$, but also the hidden data $Y_{hid}=\{x_k,k=1,\ldots,N\}$. Hence the EM is based upon a data augmentation scheme, such that the observed data are a mapping of the augmented data $Y_{obs}={\frak{m}}(Y_{aug})$, where $Y_{aug}=\{Y_{obs},Y_{hid}\}$.

The state equation (6) is non-linear, so an approximation is needed during the E-step. The likelihood and the conditional likelihood are approximated based on a linearization, and the UF is used for filtering and smoothing. The next two sub-sections briefly present the implementation of the UF and the EM algorithm for the state-space model (6)-(7).

---

4.1. The unscented filter

The UF was introduced by Julier and Uhlmann [12] as an alternative to the extended Kalman filter (EKF) for non-linear state-space models. This filter does not require the calculation of the Jacobians, and it is computationally at most as expensive as the EKF.

Both the EKF and the UF approximate the state distribution with a Gaussian one. However, instead of using the EKF linearization approach, the UF employs a deterministic sampling. The sample points completely capture the true mean and the true covariance. In contrast to the first-order accuracy of the EKF, the UF is capable to accurately capture the true posterior mean and the covariance up to the third order for a nonlinear system.

For the state-space model (6)-(7), let us denote $\bar{x}_i=E[x_i|y_1,\ldots,y_i]$ and $\hat{x}_{i+1}=E[x_{i+1}|y_1,y_2,\ldots,y_i]$ for the filtered and the predicted values, respectively, where $E[\cdot|\cdot]$ represent the conditional expectation. The corresponding conditional covariances are $\bar{P}_i=E\left [ (x_i-\bar{x}_i)(x_i-\bar{x}_i)^t|y_1,\ldots,y_i\right ]$ and $\hat{P}_{i+1}=E\biggl [ (x_{i+1}-\hat{x}_{i+1})(x_{i+1}-$ $\hat{x}_{i+1})^t|y_1,\ldots,y_i\biggl ]$. Based on the available observations $y_1\ldots y_N$, smoothed values $x_{i|N}=E[x_i|y_1$ $\ldots y_N]$ and $P_{i|N}=E[\left(x_i-x_i^N\right)\left(x_i-x_i^N\right)^t|y_1\ldots y_N]$ can be calculated.

Let denote

$f(x)=x+h\begin{bmatrix} \beta x[1](1-\frac{x[1]}{K})- \alpha x[2] x[1] \\ \lambda_1^2 x[3] -\eta_1^2x[2] \\ \lambda_2^2x[2]x[1] -\eta_2^2x[3]x[1] \end{bmatrix}, x\in {\mathbb R}^3$

(8)

and consider the linearised system

$x_{k+1}=f(\bar{x}_k)+\frac{\partial f}{\partial x}(\bar{x}_k)(x_k-\bar{x}_k)+v_{k+1},$

(9)

$y_{k}=Cx_k+w_k,.$

(10)

Notice that for the linearized system we have

$\bar{P}_{x_kx_{k+1}}=E[(x_k-\bar{x}_k)(x_{k+1}-\hat{x}_{k+1})^t|x_1,\ldots ,x_k]=\bar{P}_{k}\frac{\partial f}{\partial x}^t(\bar{x}_{k}),\label{cove}$

(11)

and with

$J_k=\bar{P}_{x_kx_{k+1}}\hat{P}_{k+1}^{-1},\label{j1}$

(12)

we have the following backward recursions

$x_{k-1|N}=\bar{x}_{k-1}+J_{k-1}(x_{k|N}-\hat{x}_k),$

(13)

$P_{k-1|N}=\bar{P}_{k-1}+J_{k-1}(P_{k|N}-\hat{P}_k)J_{k-1}^t.$

(14)

The UF filtering algorithm [11] can be used to compute $\bar{x}_{k+1}$, $\bar{P}_{k+1}$, $\hat{x}_{k+1}$, $\hat{P}_{k+1}$ and $\bar{P}_{x_kx_{k+1}}$. The smoothed values $x_{k|N}$ and $P_{k|N}$, $k=1,\ldots,N$ can be computed using (12) and the backwards recursions given in (13)-(14), starting with $k=N$.

---

4.2. The EM algorithm

The algorithm starts with an initial guess $\theta_0$ for the unknown parameters and iteratively compute the estimation $\theta^*$. Each iteration consists of two steps: the expectation (E) and the maximization (M) step. Using the current estimation $\theta_n$ of the parameters, the E-step computes the conditional expectation of the augmented data log-likelihood $Q(\theta|\theta_n)=E[\log p(\theta|Y_{aug})|$ $Y_{obs},\theta_n]$. The M-step performs a maximization with respect to the parameters $\theta$:

$\theta_{n+1}=\arg \max\limits_{\theta}Q(\theta|\theta_n).$

We briefly explain the implementation of the proposed EM-algorithm for estimating the parameters $\Theta$ of he the state-space model (6)-(7). First, we augment the data with the hidden variables $x_i$, $i=1,\ldots, N$ and we calculate the complete $\log$-likelihood:

$\log(L)=\log P(x_1,\ldots, x_N, y_1,\ldots,y_N)=\log P(y_N|y_{N-1},\ldots,y_1,x_N,\ldots, x_1)\notag\\ +\ldots+\log P(y_1|x_N,\ldots,x_1)+\log P(x_N|x_{N-1},\ldots,x_1)+\ldots+\log P(x_1).$

Using the linearized equations (9) and (10), we obtain:

$\log(L)=-5N \ln (2\pi)-\frac{1}{2}\sum\limits_{i=1}^{N}\ln(|R|)-\frac{1}{2}\sum\limits_{i=1}^{N}(y_i-Cx_i)^tR^{-1}(y_i-Cx_i)\notag\\ -\frac{1}{2}\sum\limits_{i=2}^{N}\ln (|Q|)-\frac{1}{2}\sum\limits_{i=2}^{N}(x_i-F_{i-1}x_{i-1}-b_{i-1})^tQ^{-1}(x_i-F_{i-1}x_{i-1}-b_{i-1})\notag\\ -\frac{1}{2}\ln (|\Sigma|)-\frac{1}{2}(x_1-\mu)^t\Sigma^{-1}(x_1-\mu),\notag$

where

$F_i = \begin{bmatrix} 1+h \left( \beta-2\beta \bar{x}_i[1]/K- \alpha \bar{x}_i[2] \right) &-h\alpha \bar{x}_i[1] &0\\0&1-\eta^2_1 h&\lambda^2_1 h \\h \lambda^2_2 \bar{x}_i[2]-h \eta^2_2 \bar{x}_i[3] &h\lambda^2_2\bar{x}_i[1]&1-h\eta^2_2\bar{x}_i[1] \end{bmatrix}$

and

$b_i = \begin{bmatrix} h \left( \beta\bar{x}_i[1]^{2}/K+\alpha\bar{x}_i[2]\bar{x}_i[1] \right) \\0\\ h \left( -\lambda^2_2\bar{x}_i[2]\bar{x}_i[1]+\eta^2_2\bar{x}_i[3]\bar{x}_i[1] \right) \end{bmatrix}.$

The EM algorithm iteratively maximizes $\hat{E}=E[\log(L)|y_1,\ldots,y_N]$. Using the previous formulas and $tr(AB)=tr(BA)$, where $tr$ is the notation for the trace of a matrix, we get

$\hat{E}=-5N \ln (2\pi)-\frac{N}{2}\log(|R|)-\frac{1}{2}\sum\limits_{n=1}^{N}tr[R^{-1}(y_ny_n^t-y_nx_{n|N}^tC^t-Cx_{n|N}y_n^t\notag\\ +CP_{n|N}C^t)]-\frac{1}{2}\sum\limits_{n=2}^{N}\log(|Q|)-\frac{1}{2}\sum\limits_{n=2}^{N}tr[Q^{-1}(P_{n|N}-P_{n,n-1|N}F^t_{n-1}-x_{n|N}b^t_{n-1}\notag\\ -F_{n-1}P^t_{n,n-1|N}+F_{n-1}P_{n-1|N}F^t_{n-1}+F_{n-1}x_{n|N}b^t_{n-1}-bx^t_{n|N}+b_{n-1}x^t_{n|N}F^t_{n-1}\notag\\ +b_{n-1}b^t_{n-1})]-\frac{1}{2}\log(|\Sigma|)-\frac{1}{2}tr[\Sigma^{-1}(P_1-\mu x_{1|N}^t-x_{1|N}\mu^t+\mu\mu^t)]$

The smoothed values $P_{n|N}=E[x_nx_n^t|y_1,\ldots,y_N]$, $P_{n,n-1|N}=E[x_nx_{n-1}^t|y_1,\ldots,y_N]$, and $x_{n|N}=E[x_n|$ $y_1,\ldots,y_N]$ can be calculated using the UF smoothing presented in the previous section. For the M-step, analytical update equations for the parameters $\Theta$ can be found taking derivatives with respect to the parameters $\Theta$ in the previous formula for $\hat{E}$. The implementation of the EM algorithm is summarized in Table 2

The main advantage of using the EM-algorithm is the guaranteed convergence [26,18]. However, depending on the initial guess, the algorithm can only converge to a local maximum. Various strategies for choosing the initial guess are proposed in [1], and methods for accelerating the convergence are presented in [17].

---

5. Model validation

We divide the experimental data into a training set (containing $N$ observations, where $N$ equals around 70$%$ of the data), and a test set. We use the observations in the training set to estimate the parameters. Once the parameters are estimated, we can predict the future values of $x_i$, $i=N+1,N+2, \ldots$, and we validate the model by comparing these predictions with the experimental data in the test set.

We start with the validation of the logistic model on the negative control cell growth data, and we estimate the cell growth rate $\beta$ and capacity $K$ (which we suppose are unchanged for the whole experiment). Since the logistic equation (1) can be solved analytically, the nonlinear least square method is an efficient and reliable method for parameter estimation and prediction.

Figure 4 shows the fitting of the negative control data of four replicates. The data predicted by the logistic model is in agreement with all the experimental data. Moreover the estimated parameters of all replicated experiments are consistent. Thus we conclude that the logistic equation (1) is a suitable model to represent the cell growth in non-responsive environment.

To test the responsive model (5), we present here the results obtained for PF431396, monastrol, ABT888, and HA1100 hydrochloride. These chemicals are selected from Cluster Ⅰ: DNA/RNA-nucleic acid targets, and Cluster Ⅹ: protein -motor targets. These two clusters are the top 2 most common clusters in the 63 compounds experiments, covering more than half of the cases.

We perform preprocessing of the experiment data before testing the toxicity model. Since the cell index is not measured uniformly on time, we apply a smoothing spline to approximate the experimental data. Figure 5 presents a good agreement between the experimental data and the fitting spline.

After interpolation with a smoothing spline, we sample uniformly from the approximating spline, and we use the first 70% experimental data (i.e. the training set) to estimate the parameters, $\alpha$, $\eta_1$, $\eta_2$, $\lambda_1$ and $\lambda_2$ with the EM algorithm. The estimated values of the parameters for the four chemicals mentioned before are given in Table 1. Using the estimated values of the parameters, we apply the UF to predict the cell growth/death, and we compare the predicted values with the data in the test set (i.e. with the last 30% of the experimental data).

From Figs. 6-8, the predictions are in agreement with the experimental data. We have plotted only the TCRCs corresponding to the largest 8 out of 11 values of the initial concentration of the chemicals because the curves corresponding to the first 3 values are similar to the negative control. The observations from the initial time $t_0$ to the red vertical line represent the training set, and are used for parameter estimation. The test set is formed with the data after the vertical red line. In all the figures the experimental data is plotted using dotted lines, and the filtered or predicted observations are presented using plain lines. The concentration reduces by a third from top to bottom, and from left to right.

---

6. Parameter uncertainty and sensitivity analysis

In addition to the point estimates obtained with the EM algorithm, we analyze the effect of parameter uncertainty on predicting the cell index using Monte Carlo simulations and global sensitivity analysis (GSA). Sensitivity Analysis (SA) refers to a broad group of methods that ranks parameters by their effect on output variables. GSA is the most suitable group of sensitivity analysis methods for application in highly nonlinear dynamic models [13]. In GSA, variations of all input parameters simultaneously are considered. As a result, interactions among different inputs can be detected. There is a variety of GSA methods, and here we apply the Sobol's method [23].

---

6.1. Simulations with random parameters

To account for parameter uncertainty, we run simulations assuming that all parameters are random variables with triangle distributions on intervals centered around the estimated values of the parameters given in Table 1. Based on these simulations we determine the expected cell index and the extinction probability. We consider intervals corresponding to perturbations of 10%, we run the model for 1000 realizations and we simulate for t=60 hours (similarly with the experimental settings). In Figs. 10(a), 11(a) we display some TCRCs recorded after experiments for PF431396 and ABT888. We compare these curves with the simulations results obtained using the proposed mathematical model.

In Figs. 10(b), 11(b) we plot the expected cell index and the probability of extinction against the concentration of PF431396 and ABT888. Notice that the values of the expected cell index are in good agreement with the experimental data represented in the TCRCs. The probability of extinction increases very steeply from 0 to 1 when the concentrations increase and pass a certain threshold. This illustrates the conclusions of the stability analysis presented in section 3.

From figures 6, 8, 9 we can notice that the TCRCs for PF431396, ABT888 and HA1100 hydrochloride cover both cases corresponding to persistence and to extinction. However, from figures 7 we can notice that the experiments for monastrol are done with too small initial concentrations because all TCRCs show persistence. The results of this type of simulations can be used to determine an appropriate range for the initial concentration of chemicals $CE(0)$ used in the experiments such that both values smaller and larger than the threshold between extinction and persistence are included. They extend the stability analysis done in section 3 by considering also the effects of parameter uncertainty.

---

6.2. Global sensitivity analysis

Sobol's sensitivity measures [23] that utilize the analysis of variance (ANOVA) of the model output are among the most widely used GSA methods [13] [9]. In Sobol's method, the variance of the model output is decomposed into components that result from individual parameters as well as from parameter interactions.

Consider a mathematical model represented by the square integrable function $f(\mathbf{x})$, where $\mathbf{x} = (x_1,\cdots,x_d)$ are the parameters of the model. Without loss of generality, we suppose $x_1,\cdots,x_d$ are independent, uniformly distributed random variables on $[0,1]^d$, and $f(\mathbf{x})$ has the ANOVA type of decomposition

$f(\mathbf{x}) = f_\varnothing + \sum\limits_{k} f_{[k]}(x_k)+\sum\limits_{k_1<k_2}f_{[k_1,k_2]}(x_{k_1},x_{k_2}) +\cdots+ f_{[1,\cdots,d]}(x_{1},\ldots,x_{d}) \label{gsa}$

(15)

where $f_\varnothing$ is a constant and $f_{[k_1,\ldots,k_r]}(x_{k_1},\cdots, x_{k_r})$ is a function only depending on the variables $x_{k_1},\cdots x_{k_r}$ such that for $1\leq i\leq r$,

$\int_0^1 f_{[k_1,\ldots,k_r]}(x_{k_1},\cdots, x_{k_r}) dx_{k_i} = 0.$

Thus the summands in (15) are orthogonal to each other and we have

$\int_{[0,1]^d} f(\mathbf{x}) d\mathbf{x} = f_\varnothing.$

Let us define

$D:=\int_{[0,1]^d} f^2(\mathbf{x}) d\mathbf{x} - f^2_\varnothing$

$D_{k_1,\ldots,k_r} := \int_{[0,1]^r} f^2_{[k_1,\ldots,k_r]}(x_{k_1},\cdots, x_{k_r})dx_{k_1} \cdots dx_{k_r}.$

We have Var$[f(\mathbf{x})] = D$, Var$[E[f(\mathbf{x})|x_i]] = D_i$, Var$[E[f(\mathbf{x})|x_{i},x_j]] = D_{i}+D_j+D_{i,j}$, $\cdots$, so $D$ can be seen as the total variance and $D_{k_1,\ldots,k_r}$ as partial variances. The Sobol's indices are defined as

$S_{k_1,\ldots,k_r} = \frac{D_{k_1,\ldots,k_r}}{D}.$

Squaring (15), integrating over $[0,1]^d$ and using the orthogonality of the summands [23], we get

$\sum\limits_{r=1}^d \sum\limits_{k_1<k_2<\cdots<k_r} D_{k_1,\ldots,k_r} = D.$

The GSA indices of a model can be evaluated using Monte Carlo numerical integration. To demonstrate the process, we present the estimation of first order indices $D_i$ for $1\leq i \leq d$, and the second order indices $D_{i,j}$ for $1\leq i<j\leq d$.

Let $\mathbf{\xi}_i= (\xi^i_1,\ldots,\xi^i_d)$ for $i = 1,\ldots,N$ be random samples of $\mathbf{x}$. Then

$f_\varnothing \approx \frac{1}{N}\sum\limits_{i=1}^N f(\mathbf{\xi}_i), D \approx \frac{1}{N}\sum\limits_{i=1}^N f^2(\mathbf{\xi}_i) - f^2_\varnothing .$

Then we generate another group of samples $\hat{\mathbf{\xi}}_i= (\hat{\xi}^i_1,\ldots ,\hat{\xi}^i_d)$ for $i = 1,\ldots,N$, which are independent with $\mathbf{\xi}_i$. The GSA indices are estimated as

$\begin{split} D_{j} &\approx \frac{1}{N}\sum\limits_{i=1}^N f(\xi^i_1,\ldots ,\xi^i_d)f(\hat{\xi}^i_1,\ldots,\hat{\xi}^i_{j-1},\xi^i_j,\hat{\xi}^i_{j+1},\ldots,\hat{\xi}^i_d) -f^2_\varnothing \\ D_{j_1,j_2} &\approx \frac{1}{N}\sum\limits_{i=1}^N f(\xi^i_1,\ldots ,\xi^i_d)f(\hat{\xi}^i_1,\ldots,\hat{\xi}^i_{j_1-1},\xi^i_{j_1}\hat{\xi}^i_{j_1+1},\ldots ,\hat{\xi}^i_{j_2-1},\xi^i_{j_2},\hat{\xi}^i_{j_2+1}, \ldots ,\hat{\xi}^i_d) -f^2_\varnothing \end{split}$

We present the GSA studies for our proposed model. The estimated values of the parameters for different chemicals are given in Table 3. We consider that all parameters are random variables with uniform distribution on intervals corresponding to perturbations of 10%. To calculate the Sobol's indices we use the Matlab toolbox named GSAT developed by Cannavó [2].

In Figs. (12), (13), we present the ranking for the GSA indices for PF431396 and ABT888 with various concentration, respectively. We can notice that when the chemicals' concentration is small the capacity $K$ is the most sensitive parameter of the model. This can be explained by the fact that for low concentrations of toxicant the cell population is persistent, and asymptotically the cell index $n$ approaches the capacity $K$ (see section 3). As the concentration increase to the transition stage from persistence to extinction, the cell population become more sensitive to the intake and outtake rates $\lambda_1,\eta_1,\lambda_2,\eta_2$.

Figure (14) present the second order GSA indices for the same chemicals at the transition stage from persistence to extinction. The thickness of lines of any pair of parameters represents the rank of their co-sensitivity (the thicker the line, the greater the co-sensitivity). These pictures also show that the intakes and outtakes rates are the most sensitive parameters for the cells' population model. Similar results are obtained for monastrol and HA1100 hydrochloride.

---

7. Discussion and summary

The novelty of this paper is the application of a dual methodology, including an efficient statistical method for parameter estimation and prediction of the cell index and chemicals' concentration, followed by a numerical and analytical stability study for the dynamical system modeling the cells' growth/death.

We propose a mathematical model expressed by a three dimensional system of differential equations that can be expressed in a state-space form, and we design an EM algorithm based on the unscented filter for parameter estimation and prediction. Since it requires less tunning, the EM algorithm is a more efficient alternative for parameter estimation than the MHE or the EKF [8]. The results presented here for four chemicals with different MOAs illustrate the efficiency of the proposed algorithm and validate the mathematical model. Our approach is very efficient because the model parameters can be estimated based on TCRCs provided by only two experiments, namely the negative control and an experiment with a relatively large initial concentration of the chemical. Using the mathematical model, accurate TCRCs for other concentrations can be obtained numerically.

Based on the TCRCs, clustering methods or machine learning algorithms can be applied to classify the chemicals. However, some chemicals were tested with either too low concentrations (so no significant effect on the cells growth/death is noticed) or too high concentrations (so the only effect that can be noticed is the rapid cells' death). Consequently, the TCRCs generated after this kind of experiments do not have the features necessary for statistical classification. The main practical motivation of this study was to find a way to determine the lowest concentration of toxin that can kill the cells (i.e. the value of chemical's initial concentration that separate the persistent cells' populations from those that go to extinction).

Once the parameters are estimated, the stability analysis and numerical simulations can provide a range of values for the chemical's initial concentration that give TCRCs with useful features. The asymptotic analysis of the proposed model shows that persistence occurs if the initial concentration of chemical is smaller than a certain threshold. For the persistent trajectories, the cell index tends toward the nominal carrying capacity of the population, and the internal and external chemical concentrations tend toward zero. Extinction depends on the initial concentration of toxicant and the entire past history of the population. The extinction behavior is more complex because when the cell index tends to zero, the internal and external toxicant concentrations can approach an infinite number of limiting values.

We illustrate this analytical and numerical analysis for monastrol, but the proposed model based on first order kinetics gives accurate predictions for the TCRCs corresponding to the 63 chemicals considered in this study. However, there might be limitations for different experimental settings or chemicals with different MOA.

We plan to further extend this study using uncertainty quantification methods. Although the proposed model was validated by experimental data, parameter uncertainty is inherent. To account for parameter uncertainty and calculate the probability of extinction, we have considered random perturbation of parameters' values, we have done Monte Carlo simulations, and we have calculated the Sobol's GSA indices to find out which parameters are important and sensitive to change. Based on the results of the GSA we plan to introduce a model with random parameters, for example to consider a system of stochastic differential equations similar with the ones in [15,16].

---

Acknowledgments

This work was supported by a Mathematics of Information Technology and Complex Systems (MITACS) Accelerate Clusters grant and by a grant from the Natural Sciences and Engineering Research Council of Canada.

---

Appendix A. Center manifold

Let $\begin{bmatrix} 0, CE^*\frac{\lambda_1^2}{\eta_1^2},CE^* \end{bmatrix}$ be any extinction equilibrium with $CE^*\ne \frac{\beta\eta_1^2}{\alpha\lambda_1^2}$. To translate this fixed point to the origin, we consider two new variables $y_1=C_0-CE^*\frac{\lambda_1^2}{\eta_1^2}$ and $z_1=CE-CE^*$. The corresponding model in terms of $[n, y_1, z_1]$ is given by the following system of equations

$\begin{split} \frac{dn(t)}{dt} &= \left(\beta-\frac{\alpha\lambda_1^2CE^*}{\eta_1^2}\right) n(t)-\frac{\beta n(t)^2}{K}- \alpha y_1(t) n(t) \\ \frac{dy_1(t)}{dt} &= \lambda_1^2 z_1(t) -\eta_1^2y_1(t)\\ \frac{dz_1(t)}{dt} &=\frac{CE^*(-\eta_1^2\eta_2^2+\lambda_1^2\lambda_2^2)n(t)}{\eta_1^2}+ \lambda_2^2 y_1(t)n(t) -\eta_2^2z_1(t)n(t). \end{split}\label{modelo1}$

(16)

The diagonal form of this system is

$\begin{split} \frac{dx}{dt}&=-\frac{(\alpha CE^*\lambda_1^2-\beta\eta_1^2)}{\eta_1^2}x+F_1(x,y,z)\\ \frac{dy}{dt}&=-\eta_1^2y+F_2(x,y,z)\\ \frac{dz}{dt}&=F_3(x,y,z) \end{split}\label{modelo2}$

(17)

where the new variables are

$x=\frac{CE^*(\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2)}{(\alpha CE^*\lambda_1^2-\beta\eta_1^2)}n\notag\\ y=\frac{(\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2)CE^*\lambda_1^2}{(\alpha CE^*\lambda_1^2-\eta_1^4-\beta\eta_1^2)\eta_1^2}n+y_1-\frac{\lambda_1^2}{\eta_1^2}z_1\notag\\ z=\frac{CE^*(\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2)}{(\alpha CE^*\lambda_1^2-\beta\eta_1^2)}n+z_1,$

and the non-linear terms are of the form $F_i(x,y,z)=c_1x^2+c_2xy+c_3xz$, for some coefficients $c_i$, $i=1,\ldots,3$ that can be easily determined using MAPLE. The local center manifold at $[0,0,0]$ (see chapter 2.12 in [21]) is of the form

$\begin{split} x&=h(z)=a_1z^2+a_2z^3+\ldots\\ y&=g(z)=b_1z^2+b_2z^3+\ldots. \end{split}\label{mani}$

(18)

To determine the coefficients $a_i$, $b_i$, $i=1,\ldots$, we substitute the previous expansions into the following equation (see Theorem 1, chapter 2.12 in [21])

$F_3(h(z),g(z),z)D\begin{bmatrix}h(z)\\g(z)\end{bmatrix}-\begin{bmatrix}-\frac{(\alpha CE^*\lambda_1^2-\beta\eta_1^2)}{\eta_1^2}h(z)+F_1(h(z),g(z),z)\\-\eta_1^2g(z)+F_2(h(z),g(z),z)\end{bmatrix}=0$

(19)

Simple, but long calculations give $a_i=b_i=0$, $i=1,\ldots$, so the initial system is topologically conjugate (see Theorem 2, chapter 2.12 in [21]) to the system

$\begin{split} \frac{dx}{dt}&=-\frac{(\alpha CE^*\lambda_1^2-\beta\eta_1^2)}{\eta_1^2}x\\ \frac{dy}{dt}&=-\eta_1^2y\\ \frac{dz}{dt}&=0. \end{split}\label{modelo3}$

(20)

and, in terms of the initial variables, the local center manifold is the line $\{(n, C_0, CE): n=0, \lambda_1^2 CE -\eta_1^2C_0=0\}$.

---

Appendix B. The proofs

Before we give the proofs for the lemma and theorem included in section 3, we notice that the we have the following formulas for the solutions of system (5) (see also [4]). Since equations (3) and (4) are linear we have

$C_0(t)=\lambda_1^2 e^{-\eta_1^2 t}\int_0^t CE(s) e^{\eta_1^2 s} ds\label{form1}\\ CE(t)=CE(0)\exp{\left(-\eta_2 ^2\int_0^t n(s)ds\right)}+\lambda_2^2\exp{\left(-\eta_2 ^2\int_0^t n(s)ds\right)}$

(21)

$\int_0^t C_0(s) n(s)\exp{\left(\eta_2 ^2\int_0^s n(l)dl\right )}ds, t\ge 0.$

(22)

Equation (2) is a Bernoulli equation, so we get

$n(t)=\frac{n(0)\exp{\left(\beta t-\alpha\int_0^t C_0(s)ds\right)}}{1+\frac{n(0)\beta}{K}\int_0^t\exp{\left(\beta s-\alpha\int_0^s C_0(l)dl\right)}ds}, t\ge 0.\label{form3}$

(23)

---

B.1. Proof of lemma 3.1

Assume that there exists $T_1>0$ such that $CE(T_1)>CE(0)$ or $CE(T_1)\le 0$, and suppose that $[0,T]$ is the maximum interval for which $0<CE(t)\le CE(0)$. Then we should either have $CE(T)=CE(0)$ or $CE(T)=0$.

From (21) we get for any $0<t\le T$

$0\le C_0(t)=\lambda_1^2 e^{-\eta_2^2 t}\int_0^t CE(s) e^{\eta_1^2 s} ds\le \lambda_1^2 e^{-\eta_2^2 t}CE(0)\int_0^t e^{\eta_1^2 s} ds\notag\\ =\frac{\lambda_1^2 CE(0)}{\eta_1^2}\left(1-e ^{-\eta_1^2t}\right)\le \frac{\lambda_1^2 CE(0)}{\eta_1^2}$

Thus, for any $0<t\le T$

$0\le \int_0^t C_0(s) ds\le \frac{\lambda_1^2 CE(0)t}{\eta_1^2}<\infty\label{intor}$

(24)

and from equation (23) we have $0<n(t)\le n(0)e^{\beta t}$. Using this and equation (22) we get for any $0<t\le T$

$CE(t)\ge CE(0)\exp{\left(-\eta_2 ^2\int_0^t n(s)ds\right)}>0,$

(25)

so $CE(T)>0$.

Moreover, if $CE(T)=CE(0)$, then from equation (4) we obtain

$\frac{dCE}{dt}\biggl |_{t=T}=\lambda_2^2C_0(T)n(T)-\eta_2^2 CE(T)n(T)\notag\\ \le CE(0) n(T)\left(\frac{\lambda_1^2 \lambda_2^2}{\eta_1^2}-\eta_2^2\right)<0$

(26)

So there exists $\delta>0$ such that $CE(t)\le CE(0)$ for all $t\in [T,T+\delta]$, and we have a contradiction with the definition of the interval $[0,T]$.

Thus, we have proved that $0<CE(t)\le CE(0)$ for any $t\ge 0$. Using this and (21) we also get $0\le C_0(t)\le \frac{\lambda_1^2 CE(0)}{\eta_1^2}$, for any $t\ge 0$. Consequently inequalities (24) are true for any $t\ge 0$, and from (23) we get $n(t)>0$, for all $t\ge 0$.

---

B.2. Proof of Theorem 3.2

Since from Lemma 3.1 we know that $n(t)>0$, for any $t\ge 0$, we have $\int_0^\infty n(t)dt=\infty$ or $0\le \int_0^\infty n(t)dt<\infty$.

If $|n|_1=\int_0^\infty n(t)dt<\infty$ then for any $t\ge 0$ we also have

$\int_0^t C_0(s) n(s)\exp{\left(\eta_2 ^2\int_0^s n(l)dl\right )}ds\le \frac{\lambda_1^2 CE(0)}{\eta_1^2}\exp{\left(\eta_2 ^2|n|_1\right )}|n|_1,$

(27)

so $0\le\int_0^\infty C_0(s) n(s)\exp{\left(\eta_2 ^2\int_0^s n(l)dl\right )}ds=M<\infty$. Thus from (22) we get

$\lim\limits_{t\rightarrow\infty} CE(t)=CE(0)\exp(-\eta_2^2|n|_1)+\lambda_2^2 M\exp(-\eta_2^2|n|_1)<\infty.$

(28)

Consequently, there exists $T_1>0$ such that tor any $t>T_1$ we have $CE(t)>CE(0)\exp(-\eta_2^2|n|_1)/2$. This implies that $\int_0^\infty CE(s) e^{\eta_1^2s}ds=\infty$ because for any $t> T_1$ we have

$\int_0^t CE(s) e^{\eta_1^2s}ds\ge \int_{T_1}^t CE(s) e^{\eta_1^2s}ds\ge CE(0)\exp(-\eta_2^2|n|_1)/2\int_{T_1}^t e^{\eta_1^2s}ds.$

So we can apply L'Hospital rule in (21), and we get

$\lim\limits_{t\rightarrow\infty} C_0(t)=\frac{\lambda_1^2}{\eta_1^2}\lim\limits_{t\rightarrow\infty} CE(t).$

(29)

Thus, if $\int_0^\infty n(t)dt<\infty$, then $\lim_{t\rightarrow\infty} CE(t)$ and $\lim_{t\rightarrow\infty} C_0(t)$ exist and they are related by the previous equation.

Next we consider the case when $\int_0^\infty n(t)dt=\infty$. If $0\le\int_0^\infty C_0(s) n(s)$ $\exp \bigl(\eta_2 ^2$ $\int_0^s n(l)dl\bigl)ds<\infty$, from (22) we get $\lim_{t\rightarrow\infty} CE(t)=0$. On the other hand, if $\int_0^\infty C_0(s) n(s)\exp{\left(\eta_2 ^2\int_0^s n(l)dl\right )}ds=\infty$, from L'Hospital rule in (22) we have

$0\le \frac{\lambda_2^2}{\eta_2^2}\liminf\limits_{t\rightarrow\infty}C_0(t)\le \liminf\limits_{t\rightarrow\infty}CE(t)\le \limsup\limits_{t\rightarrow\infty}CE(t)\le \frac{\lambda_2^2}{\eta_2^2}\limsup\limits_{t\rightarrow\infty}C_0(t)$

(30)

Similarly, from (21) we either get that $\lim_{t\rightarrow\infty} C_0(t)=0$ (if $\int_0^\infty CE(s) e^{\eta_1^2s}ds<\infty$), or we have

$0\le \frac{\lambda_1^2}{\eta_1^2}\liminf\limits_{t\rightarrow\infty}CE(t)\le \liminf\limits_{t\rightarrow\infty}C_0(t)\le \limsup\limits_{t\rightarrow\infty}C_0(t)\le \frac{\lambda_1^2}{\eta_1^2}\limsup\limits_{t\rightarrow\infty}CE(t),$

(31)

(if $\int_0^\infty CE(s) e^{\eta_1^2s}ds=\infty$). All these possible cases give

$\lim\limits_{t\rightarrow\infty}C_0(t)=\lim\limits_{t\rightarrow\infty}CE(t)=0,$

(32)

because $\eta_1^2\eta_2^2-\lambda_1^2\lambda_2^2>0$. Thus, if $\int_0^\infty n(t)dt=\infty$, then $\lim_{t\rightarrow\infty} CE(t)$ and $\lim_{t\rightarrow\infty} C_0(t)$ exist and they are equal with zero.

In conclusion, we have shown that in both cases $\lim_{t\rightarrow\infty} CE(t)$ and $\lim_{t\rightarrow\infty} C_0(t)$ exist, and we have $\lim_{t\rightarrow\infty} C_0(t)=\frac{\lambda_1^2}{\eta_1^2}\lim_{t\rightarrow\infty} CE(t)$.

To prove the first part of the theorem we first show that if $\lim_{t\rightarrow\infty} CE(t)<\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$ then $\int_0^\infty n(t)dt=\infty$. From the previous discussion we know that $\lim_{t\rightarrow\infty} CE(t)<\frac{\beta\eta_1^2}{\alpha\lambda_1^2}$ implies $\lim_{t\rightarrow\infty} C_0(t)<\frac{\beta}{\alpha}$, and consequently

$\lim\limits_{t\rightarrow\infty}\exp{\left(\beta t-\alpha\int_0^t C_0(s)ds\right)}=\infty, \lim\limits_{t\rightarrow\infty}\int_0^t\exp{\left(\beta s-\alpha\int_0^s C_0(l)dl\right)}ds=\infty$

Hence, from L'Hospital rule in (23) we have

$0\le \lim\limits_{t\rightarrow\infty}n(t)=\frac{K\left(\beta-\alpha \lim_{t\rightarrow\infty}C_0(t)\right) }{\beta}\label{eqn0}$

(33)

If $\int_0^\infty n(t)dt<\infty$ then we should have $\lim_{t\rightarrow\infty}n(t)=0$, and this implies $\lim_{t\rightarrow\infty} C_0(t)$ $=\frac{\beta}{\alpha}$. But $\lim_{t\rightarrow\infty} C_0(t)<\frac{\beta}{\alpha}$, so $\int_0^\infty n(t)dt=\infty$. Hence $\lim_{t\rightarrow\infty}C_0(t)=\lim_{t\rightarrow\infty}CE(t)$ $=0$, and replacing in (33) we get $\lim_{t\rightarrow\infty}n(t)=K$. This completes the proof of the first part of the theorem.

Next, if $\lim_{t\rightarrow\infty} CE(t)=CE^*>\frac{\beta\eta_1^2}{\alpha\lambda_1^2}>0$, then from the previous discussion we know that $|n|_1=\int_0^\infty n(t)dt<\infty$. Moreover, $\lim_{t\rightarrow\infty} C_0(t)=\frac{\lambda_1^2}{\eta_1^2}CE^*>\frac{\beta}{\alpha}>0$, so

$\lim\limits_{t\rightarrow\infty}\exp{\left(\beta t-\alpha\int_0^t C_0(s)ds\right)}=0.$

Thus $\lim_{t\rightarrow\infty}n(t)=0$ because from (23) we have

$0<n(t)\le n(0)\exp{\left(\beta t-\alpha\int_0^t C_0(s)ds\right)}.$

---