Patient-specific parameter estimates of glioblastoma multiforme growth dynamics from a model with explicit birth and death rates

1.   Introduction

Phytoplankton-derived carbon (internal primary production) is known to be essential for the somatic growth and reproduction of zooplankton and fish. For decades, aquatic food webs were taken as systems where carbon transfer was linear from phytoplankton to zooplankton to fish ^[1].

With the application of isotope labeling, substantial research shows that terrestrial organic matter (TOM) plays an important role in the lake food webs. TOM can not only affect the lake ecosystem in physical and chemical way but also can be exploited by consumer as a resource ^[2]. And particulate organic matter (POM) of terrestrial origin can be the key factor controlling whole-lake productivity in lakes where phytoplankton productivity is low ^[1,3]. Besides, in lakes, the inputs of TOM often equal to or exceed internal primary production ^[4].

Zooplankton composed of cladocerans, cyclopoids, and calanoids represents a critical food chain link between phytoplankton and larger fish in lakes. Lot of experimental results show that zooplankton can directly consume allochthonous POM that either entered the lake in particulate form or was formed through flocculation of allochthonous dissolved organic matter (DOM) ^[5]. Nevertheless, it is mentioned that Daphnia only supported by allochthonous carbon can survive and give birth to offspring ^[6]. For the reason that phytoplankton and zooplankton are important components of lakes, plankton mechanism considering TOM is widely studied. Lots of experiments have been done by isotope labeling ^[5,7,8].

Mathematical modelling is an important tool to investigate the plankton mechanism ^{[9,10,11,12,13,14]}. Nutrient-phytoplankton-zooplankton models were studied by many researchers ^[15,16,17]. In order to enhance the knowledge of plankton mechanism, we propose phytoplankton-zooplankton models considering TOM. We focus on the ecological function of TOM in the plankton mechanism. The results of the deterministic model can fit well with some experimental results and two hypotheses supported by experimental results in ^[18]. The related two hypotheses are given in the following part:

(Ⅰ) Catchment deposition hypothesis: Allochthony which is the portion of a zooplankter's body carbon content that is of terrestrial origin increases as more TOM is exported from the surrounding catchment.

(Ⅱ) Algal subtraction hypothesis: Allochthony increases with the availability of TOM, where algal production becomes limited by shading more than it benefits from the nutrients associated with TOM.

The rest of the paper is organized as follows: in the next section, by taking TOM into account, formulation of two deterministic models are discussed. In Section 3, global dynamical properties of the three-dimensional ODE model are completely established using Lyapunov function method. In Section $4$, by including stochastic perturbation of the white noise type, we develop a stochastic model and show that there is a stationary distribution in the stochastic model. In Section $5$, the numerical simulation is given to verify the theoretical predictions. In the biological sense, the numerical results are analyzed in detail. In the last section, the conclusion is given.

2.   Deterministic models

The research on the interaction between phytoplankton and zooplankton is important. The phytoplankton is not only the producer of the marine ecosystem but also the base of every food web. Phytoplankton is consumed by zooplankton, which is resource for consumers of higher tropical levels. By transferring organic matter and energy into higher tropical levels, zooplankton plays a critical role in lake ecosystem. Usually, the interaction between the phytoplankton and zooplankton can be described by the Prey-Predator system,

Here, $P(t)$ and $Z(t)$ are the population densities of phytoplankton and zooplankton, respectively. The phytoplankton grows with rate $r_{1}$. The interaction term taken here is Holling type-Ⅰ with $b$ and $c$ as the respective rates of grazing and biomass conversion by the zooplankton for its growth such that $b > c$. The natural mortality of zooplankton is $r_{2}$. $a_{1}$ reflects the density dependence of phytoplankton. The density dependence of predator species is wildly studied because of the environment factor ^[14]. In lakes where phytoplankton productivity is low, zooplankton will compete each other for the resources to survive, which means that it is more realistic to consider the density dependence of zooplankton population in system (2.1). Furthermore, the zooplankton is resource for higher tropical levels. The mortality caused by higher trophic levels is generally modeled by a nonlinear term ^[9]. Hence, $a_{2}$ reflects the density dependence of zooplankton and the mortality of zooplankton caused by higher tropic. Here, it is easy to know that system (2.1) has two equilibria: $O(0, 0)$, $A(\frac{r_{1}}{a_{1}}, 0).$ Besides, when $cr_{1}-a_{1}r_{2} > 0$, there is a positive equilibrium $E^{*}(P^{*}, Z^{*})$, where $P^{*} = \frac{r_{1}a_{2}+r_{2}b}{a_{1}a_{2}+bc}$, $Z^{*} = \frac{cr_{1}-a_{1}r_{2}}{a_{1}a_{2}+bc}$, which is globally asymptotically stable. If $cr_{1}-a_{1}r_{2}\leq 0$, equilibrium $A(\frac{r_{1}}{a_{1}}, 0)$ is globally asymptotically stable.

Model $(2.1)$ ignores TOM in the lake ecosystem. Since TOM has critical influence in the lake food webs, it is more realistic to incorporate it into the phytoplankton-zooplankton model. Based on the important role of TOM in lake ecosystem, we propose the terrestrial organic matter-phytoplankton-zooplankton model. The meaning of the parameters and some assumptions are proposed in following part, firstly.

$C(t)$ is the concentration of TOM. $P(t)$, $Z(t)$ are the population densities of phytoplankton and zooplankton, respectively. $s$ is the constant input rate of TOM and $\delta$ is its sinking rate. The interaction term between the TOM and zooplankton is Holling type-Ⅰ. $\beta_{0}$ is TOM uptake rate by zooplankton. TOM is converted for the growth of zooplankton with rate $\beta_{3}$. The phytoplankton grows with growth rate $r$. $K$ is the carrying capacity. The interaction term between the phytoplankton and zooplankton is also Holling type-Ⅰ with $\beta_{1}, \beta_{2}$ as the respective rates of grazing and biomass conversion by the zooplankton for their growth. $\gamma$ is the natural mortality of zooplankton. $m$ measures the strength of competition among zooplankton and the mortality of zooplankton caused by higher tropic. Above statements can be seen clearly in Figure 1.

Based on the above statements, the following ordinary differential equations can be derived,

In the biological context, the population densities of phytoplankton and zooplankton are nonnegative at any time $t$. The concentration of TOM is also nonnegative at any time $t$. We get the following initial conditions,

Proposition 2.1. The solution $(C(t), P(t), Z(t))$ of $(2.2)$ with the initial condition $(2.3)$ is nonnegative and bounded in $R_{+}^{3}$.

Proof. From the equations of (2.2), it is easy to get that $C(t)\geq C(0)e^{-\int_{0}^{t}\beta_{0}Z(s)+\delta ds}\geq 0$, $P(t) = P(0)e^{\int_{0}^{t}r\left(1-\frac{P(s)}{K}\right)-\beta_{1}Z(s)ds}\geq 0$ and $Z(t) = Z(0)e^{\int_{0}^{t}\beta_{2}P(s)+\beta_{3}C(s)-\gamma-mZ(s)ds}\geq 0$. Hence, the solutions are non-negative. Clearly, $\frac{dC}{dt}\leq -\delta C+s$ yielding $C(t)\leq \max\left \{P_{0}, \frac{s}{\delta}\right \}: = C_{max}$. $\frac{dP}{dt}\leq rP(1-\frac{P}{K})$ yielding $P(t)\leq \max\{K, P_{0}\}: = P_{max}$. Therefore, $Z(t)\leq \beta_{2}P_{max}Z+\beta_{3}C_{max}Z-mZ^{2}$. Hence, $Z(t)\leq \max \left \{Z_{0}, \frac{\beta_{2}P_{max}+\beta_{3}C_{max}}{m}\right \}.$ It implies the solutions $(C(t), P(t), Z(t))$ of $(2.2)$ are nonnegative and bounded for all $t\geq 0.$

The nonnegativity of the solution ensures that the model has practical significance. In the following part, we can get the equilibria of $(2.2)$. By direct calculation, it can be known that $(2.2)$ always has two equilibria,

Proposition 2.2. If $s > \frac{\delta \gamma}{\beta_{3}}$, there exists unique semi-trivial equilibrium $E_{3}(C_{3}, 0, Z_{3})$ of (2.2). If $K\geq \frac{\gamma}{\beta_{2}}$ and $s\in(0, s_{1})$, or $K < \frac{\gamma}{\beta_{2}}$ and $s\in(s_{0}, s_{1})$, where

there exists a unique positive equilibrium $E_{4}(C_{4}, P_{4}, Z_{4})$.

Proof. The equilibria satisfy the following system of algebraic equations,

If $P = 0$, $Z\neq 0$, the equilibrium is given by

where $C^{*}$ and $Z^{*}$ are constants. From the second equation of (2.5), we can get

where $a_{3} = -m\beta_{0}, b_{3} = -(\gamma\beta_{0}+\delta m), c_{3} = \beta_{3}s-\gamma\delta$. If $s > \frac{\delta \gamma}{\beta_{3}}$, then $c_{3} > 0$. Thus, $\Delta = b_{3}^{2}-4a_{3}c_{3} > 0$, $Z_{31}+Z_{32} = -\frac{b_{3}}{a_{3}} < 0$, $Z_{31}Z_{32} = \frac{c_{3}}{a_{3}} < 0$, which means that (2.6) exists only one positive root. It is easy to get that $Z_{3} = \frac{-b_{3}+\sqrt{b_{3}^{2}-4a_{3}c_{3}}}{2a_{3}}$. Then $C_{3} = \frac{\gamma+mZ_{3}}{\beta_{3}}.$

If $P\neq0, Z\neq0,$ the equilibria of $(2.4)$ are given by

where

If $s < s_{1}$, then $c_{4} < 0$. Hence, $H(0) = c_{4} < 0$. If $K\geq \frac{\gamma}{\beta_{2}}$, or $K < \frac{\gamma}{\beta_{2}}$ and $s > s_{0}$, $H(K) = \beta_{3}s+\beta_{2}K\delta-\gamma\delta > 0.$ And $\lim_{P \to +\infty} H(P) < 0.$ Therefore, function $H(P)$ has two positive roots, one of which is larger than $K$. However, $Z = \frac{r}{\beta_{1}}\left(1-\frac{P}{K}\right)$. The root which is bigger than $K$ should be discarded. Hence, there is a unique positive equilibrium.

Here, in biologic context, $E_{1}$ means that both the phytoplankton and zooplankton go extinct eventually. $E_{2}$ means that the zooplankton goes extinct, while the phytoplankton will exist at the density of $K$ eventually. $E_{3}$ means that the phytoplankton goes extinct, while the zooplankton will exist at the density of $Z_{3}$ eventually. $E_{4}$ means that the phytoplankton and zooplankton coexist. The existence of $E_{3}$ and $E_{4}$ is mainly determined by $s$, the constant input rate of TOM.

In the next section, we will analyse the stability of those equilibria. It should be pointed out that $(2.2)$ is totally different from the nutrient-phytoplankton-zooplankton model, for the reason that $s$ stands for TOM, which is ingested by zooplankton.

3.   Dynamical properties of (2.2)

In this section, global dynamical properties of $(2.2)$ are established.

Firstly, at the trivial equilibrium point $E_{1}(\frac{s}{\delta}, 0, 0)$ in (2.2), the corresponding characteristic equation is

It is apparent that equation (3.1) has one positive eigenvalue, which implies that $E_{1}(\frac{s}{\delta}, 0, 0)$ is always unstable.

Theorem 3.1. If $s\geq s_{1}$, then $E_{3}(C_{3}, 0, Z_{3})$ is globally asymptotically stable.

Proof. Consider the following Lyapunov function

Since the function $h(z) = z-1-\ln z$ is always nonnegative, it is easy to know that $V_{1}(t)$ is also nonnegative at any time t. Then the derivative of $V_{1}(t)$ along the solution of $(2.2)$ is given by

If $s\geq s_{1}$, then $\beta_{2}\left(\frac{r}{\beta_{1}}-Z_{3}\right)\leq 0.$ It follows that $\frac{dV_{1}}{dt}\leq 0$. Furthermore, $\frac{dV_{1}}{dt} = 0$ if and only if $C(t) = C_{3}$, $P(t) = 0$ and $Z(t) = {Z_{3}}$. Thus, the largest invariant set on which $\frac{dV_{1}}{dt}$ is zero consists of just equilibrium $E_{3}$. Therefore, by LaSalle's Invariance Principle, equilibrium $E_{3}$ is globally asymptotically stable.

Theorem 3.2. If $K\geq \frac{\gamma}{\beta_{2}}$ and $s\in(0, s_{1})$, or $K < \frac{\gamma}{\beta_{2}}$ and $s\in(s_{0}, s_{1})$, then $E_{4}(C_{4}, P_{4}, Z_{4})$ is globally asymptotically stable.

Proof. Consider the following Lyapunov function

Similarly, the derivative of $V_{2}(t)$ along the solution of $(2.2)$ is given by

It is easy to know $\frac{dV_{2}}{dt}\leq 0$. $\frac{dV_{2}}{dt} = 0$ if and only if $C(t) = C_{4}$, $P(t) = P_{4}$ and $Z(t) = Z_{4}$, and hence the largest invariant set of (2.2) in the set $\left\{ (C(t), P(t), Z(t))|\frac{dV_{2}}{dt} = 0 \right\}$ is the singleton $\{E_{4}\}$. Therefore, by the LaSalle's Invariance Principle, equilibrium $E_{4}$ is globally asymptotically stable.

Theorem 3.3. If $K < \frac{\gamma}{\beta_{2}}$ and $s\in(0, s_{0}]$, then $E_{2}(\frac{s}{\delta}, K, 0)$ is globally asymptotically stable.

Proof. Consider the following Lyapunov function

Similarly, the derivative of $V_{3}(t)$ along the solution of $(2.2)$ is given by

If $s\leq\frac{\delta \gamma}{\beta_{3}}-\frac{\beta_{2}K\delta}{\beta_{3}}$, then $\beta_{3}C_{2}+\beta_{2}P_{2}-\gamma = \frac{\beta_{3}s}{\delta}+\beta_{2}K-\gamma\leq 0$. It is easy to know $\frac{dV_{3}}{dt}\leq 0$. $\frac{dV_{3}}{dt} = 0$ if and only if $C(t) = \frac{s}{\delta}$, $P(t) = K$ and $Z(t) = 0$. Similarly, by the LaSalle's Invariance Principle, every solution of $(2.2)$ tends to $M$, where $M = \left\{ (C(t), P(t), Z(t))|\frac{dV_{3}}{dt} = 0 \right\} = \{E_{2}\}$ is the largest invariant set of (2.2). It implies that $E_{2}$ is globally asymptotically stable.

Global dynamical properties of $(2.2)$ are shown clearly in Table 1. From biological viewpoint, Theorem 3.1 implies that if the constant input of TOM is larger than $s_{1}$, it will result in the extinction of phytoplankton eventually which is impossible in the system (2.1). It also means that the increasing input of the TOM will inhibit the growth of phytoplankton, which fits well with the hypotheses(Ⅱ). In addition, in the extinct process of the phytoplankton, in order to survive, the zooplankton will uptake more TOM which is easy to get. That is the Allochthony of zooplankton will increase, which fits well with the hypotheses(Ⅰ). The coexistence of TOM and zooplankton is mentioned in ^[6]. It also means that the producers of catchment can support both the landscape and lake ecosystem, as told in ^[2]. Furthermore, from Theorem 3.2, it can be known that different from the system $(2.1)$, with the enough constant input of the TOM($s > s_{0}$), zooplankton in the deterministic setting can produce and give offsprings at small carrying capacity(K). The result fits well with the experimental result in ^[19]. Besides, the threshold value $s_{1}$ is positively related to the birth rate of phytoplankton, which means that phytoplankton with high birth rate is less possible to become extinct. The result of Theorem 3.3 implies that zooplankton goes extinct if both the carrying capacity and constant input rate of TOM are small.

4.   Mathematical model with stochastic perturbation

In lake ecosystem, the weather, temperature and other physical factors are hard to predict. TOM, phytoplankton and zooplankton are easily affected by these factors. The deterministic model studied in the previous section does not take environmental randomness into consideration. Taking the randomness into account is more realistic. Recently, many researchers studied the stochastic biological models. Different types of stochastic perturbations are introduced into the models. Prey-predator models with fluctuations around the positive equilibrium were studied ^[20,21]. The stochastic perturbations which are proportional to the variables were discussed ^[22,23]. Fluctuations manifesting in the transmission coefficient rate were discussed ^[24,25].

In this section, by considering role of unpredictable environmental factors in plankton dynamics, we get following stochastic equations,

where $B_{i}(t), i = 1, 2, 3$, are independent Brownian motions and $\sigma_{i}$, $i = 1, 2, 3$ are the corresponding intensities of stochastic perturbations. Obviously, all the equilibria of $(2.2)$ are no longer the equilibria of the stochastic model (4.1).

In this paper, let $(\Omega, \mathcal{F}, \{\mathcal{F}_{t}\}_{t\geq 0}, \mathrm{P})$ be a complete probability space with filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ satisfying the usual conditions(i.e. it is right continuous and $\mathcal{F}_{0}$ contains all $\mathrm{P}$-null sets). Let $B_{i}(t)$, $i = 1, 2, 3,$ be a standard one-dimensional Brownian motion defined on this complete probability space. Let $R^{3}_{+} = \{x \in R^{3}: x_{i} > 0, 1\leq i \leq 3\}$.

Definition 4.1 ([26]). A non-negative random variable $\tau(\omega)$, which is allowed to the value $\infty$, is called stopping time(with respect to filtration $\mathcal{F}_{t}$), if for each t, the event $\{\omega: \tau \leq t\}\in \mathcal{F}_{t}$.

Lemma 4.2 ([26]). Let $\{X_{t}\}_{t}$ be right-continuous on $R^{n}$ and adapted to $\mathcal{F}_{t}$. $Define \inf \emptyset = \infty$. If D is an open interval on $R^{n}$, $\tau_{D} = \inf\{t\geq 0: X_{t}\not\in D\}$ is a $\mathcal{F}_{t}$ stopping time.

Lemma 4.3 ([26]). Let $x(t) = (x_{1}(t), ..., x_{n}(t))$ be a regular adapted process $n$-dimensional vector process. We can have any number $n$ of process driven by a $d$-dimensional Brownian motion

where $\sigma(t)$ is $n\times d$ matrix valued function, $B_{t}$ is $d$-dimensional Brownian motion, $x(t), \; b(t)$ are $n$-dim vector-valued functions, the integrals with respect to Brownian motion are It$\hat{o}$ integrals. Then $x(t)$ is called an It$\hat{o}$ process. The only restriction is that this dependence results in:

(i) for any $i = 1, 2...n$, $b_{i}(t)$ is adapted and $\int_{0}^{T}\left|b_{i}(t)\right | < \infty$ a.s..

(ii) for any $i = 1, 2, ...n$, $\sigma_{ij}(t)$ is adapted and $\int_{0}^{T}\sigma_{ij}^{2}(t) < \infty$ a.s.

Lemma 4.4 ([27]). Let $x(t)$ be a $l$-dimensional $It\hat{o}$ process on $t\geq 0$ with stochastic differential

where $f \in \mathcal{L}(R_{+}; R^{l})$ and $g \in \mathcal{L}^{2}(R_{+}; R^{l\times m}).$ Let $V\in C^{2, 1}(R^{l}\times R_{+}; R)$. Then $V(x(t), t)$ is till an $It\hat{o}$ process with the stochastic differential given by

where

Suppose $V\in C^{2, 1}(R^{l}\times R_{+}; R)$, define a operator $LV$ from $R^l\times R_{+}$ to $R$ by

Theorem 4.5. For any given initial value $(C(0), P(0), Z(0))\in R^{3}_{+}$, (4.1) has almost surely (a.s.) a unique positive solution $(C(t), P(t), Z(t))$ for $t\geqslant 0$, and the solution remains in $R^{3}_{+}$ with probability 1.

Proof. Since the coefficients of the equation are locally Lipschitz continuous, for any given initial value $(C(0), P(0), Z(0))\in R^{3}_{+}$, it can be obtained that there is a unique local solution $(C(t), P(t), Z(t))$ on $t\in[0, \tau_e)$ where $\tau_{e}$ is the explosion time. Explosion refers to the situation when the process reaches infinite values in finite time. Solution can be considered until the time of explosion^[26]. If we can check that $\tau_e = \infty$ a.s., then the solution is global. Let $h_{0}$ be sufficient large, such that every coordinate of (C(0), P(0), Z(0)) lies within the interval $\left[\frac{1}{h_{0}}, h_{0}\right]$. For each integer $h\geq h_{0}$, stopping time is defined by

where we set $\inf \emptyset = \infty.$ Clearly, $\tau_{h}$ is increasing as $h\rightarrow\infty$. Let $\tau_{\infty} = \lim_{h\rightarrow\infty}\tau_{h}$, where $\tau_{\infty} \leq \tau_{e}$ a.s. If $\tau_{\infty} = \infty$ a.s., then $\tau_{e} = \infty$ and $(C(t), P(t), Z(t))\in R^{3}_{+}$ a.s. for all $t\geq 0$. If not, then it exists constants $T > 0$ and $\epsilon\in(0, 1)$ such that $P\{\tau_{\infty}\leq T\} > \epsilon$. Hence, there is integer $h_{1}\geq h_{0}$ such that

Define a function $V_{4}$: $R^{3}_{+}\rightarrow R_{+}$ as

We can get the non-negativity of this function from $f(u) = u+1-\ln u \geq 0, \forall u > 0$. Applying It$\hat{o}$ formula, it follows that

Hence,

Let

Since the function $f$ is positive, it yields $u_{i}\leq 2(u_{i}+1-\ln{u_{i}})$^[23]. Therefore, $\beta_{3}C+\frac{\beta_{2}}{\beta_{1}}\left(r+\frac{r}{K}+\beta_{2}\right)P+q_{3}Z\leq 2q_{2}V_{4}(C, P, Z)$, which implies

and hence,

where $q_{4} = \max\{q_{1}, 2q_{2}\}$. Let us define $a\wedge b = \min\{a, b\}.$ If $h_{1}\leq T$, then

Taking expectation of both sides, we can get that

By the Gronwall inequality,

We set $\Omega_{h} = \{\tau_{h}\leq T\}$ for $h\geq h_{1}$, which implies $P(\Omega_{h})\geq\epsilon$. Note that for every $\omega\in\Omega_{h}$, we can have that $C(\tau_{h}, \omega)$ equals either $h$ or $\frac{1}{h}$ or $P(\tau_{h}, \omega)$ equals either $h$ or $\frac{1}{h}$ or $Z(\tau_{h}, \omega)$ equals either $h$ or $\frac{1}{h}$, and hence,

Then it implies that

where $1_{\Omega_{h}}$ is the indicator function of $\Omega_{h}$. Then, we let $h \to \infty$, which leads to the contradiction $\infty > (V_{4}(C(0), P(0), Z(0))+q_{4}T)e^{q_{4}T} = \infty$. Therefore, $\tau_{h} = \infty$ a.s.

Now, we give the main result of this section in the following part.

Theorem 4.6. Assume $\sigma_{1}^{2} < \delta$, $\sigma_{2}$, $\sigma_{3} > 0$ such that

Then there is a stationary distribution $\mu(\cdot)$ for (4.1) with initial value $(P_{0}, C_{0}, Z_{0})\in R_{+}^{3}$, which has ergodic property.

Proof. Define a Lyapunov functional

Then applying It$\hat{o}$'s formula to (4.1), we obtain

Since $(a+b)^{2}\leq 2(a^{2}+b^{2})$, we have that

where $\eta = \frac{\beta_{3}\sigma_{1}^{2}}{\beta_{0}}C_{4}+\frac{\beta_{2}\sigma_{2}^{2}}{2\beta_{1}}+\frac{\sigma_{3}^{2}}{2}$. When $\eta < \min \left\{\frac{\beta_{3}}{\beta_{0}}(\delta-\sigma_{1}^{2})C_{4}, \; mZ_{4}^{2}, \; \frac{r\beta_{2}}{K\beta_{1}}P_{4}^{2}\right\}$, we can get that the ellipsoid

lies entirely in $R_{+}^{3}$. Let $U$ be a neighborhood of the ellipsoid with $\bar{U}\subset E_{l} = R_{+}^{3}$, so for $x \in U\backslash E_{l}$, $LV_{5} < -\zeta$($\zeta$ is a positive constant). Therefore, we have that condition (B.2) in Lemma $2.1$ of ^[23] is satisfied($E_{l}$ denotes euclidean $l$-space). Besides, there is $M > 0$ such that

Thus, we have that the condition (B.1) of ^[23] is also satisfied. Therefore, there exists a stable stationary distribution $\mu(\cdot)$ which is ergodic in $(4.1)$.

5.   Numerical simulations

In this section, the numerical simulations of $(2.2)$ and $(4.1)$ are given to study the proposed models in details. Different simulated results can fit with different hypotheses mentioned in the first section. Unless otherwise stated, parameter values given in Table 2 are used for the simulations. In order to get biologically plausible results, many values are taken from parameter ranges found in the literature. And the numerical simulations of (2.2) are prepared by tool kit ode45 of Matlab.

It is easy to get $s_{1} = 4.4,$ which is the threshold value of model $(2.2)$. It should be noticed that since it is impossible that the absorbed energy is fully used for reproduction, the values of our parameters should satisfy $\beta_{0} > \beta_{3}$ and $\beta_{1} > \beta_{2}$. In the experimental results of ^[13], it points out that for zooplankton, the growth efficiency of assimilating TOM is lower than that of ingesting phytoplankton. In following analysis, hence, the values of our parameters should satisfy $\beta_{1} > \beta_{0}.$

In the case where $s = 1 < s_{1}$, Figure 2(a) shows that solutions dampen and tend to a stable steady state. The carrying capacity varies with the geographical location of the lake. Especially, the carrying capacity of poor nutrient lake is very small. If we decrease $K$ to 1, from Theorem $3.1,$ we can know dynamical behavior of (2.2) is similar. However, the solution of (2.1) with the same parameter values will tend to $A(\frac{r_{1}}{a_{1}}, 0)$, which means the extinction of zooplankton. Above statements can fit well with the experimental results of ^[28] that Daphnia magna(zooplankton) can use terrestrial-derived dissolved organic matter (t-DOM) to support growth and reproduction when alternative food sources are limiting.

The catchment areas around different lakes are different. The constant input of TOM is mainly determined by the catchment areas. It is reasonable to consider the variation of parameter value $s.$ The values of other parameters are kept invariant.

In the case where $s = 5 > s_{1}$, we see from Figure 3(a) that the phytoplankton goes extinct eventually. In particular, because of the extinction of phytoplankton, the growth of zooplankton is determined by TOM. That is, the zooplankton shows a preference for TOM, which fits well with the hypothesis (Ⅰ).

In addition, it also means the increasing input rate of TOM will hinder the growth of phytoplankton, which fits well with the hypothesis (Ⅱ).

In this part, we decrease $K = 100$ to $1$ and $s = 1$ to $0.1$. The other parameter values are the same as before. The result of numerical simulation is shown in Figure 4(a). We can notice that the density of zooplankton decreases to zero, which means if both environmental capacity and constant input of TOM are small, the zooplankton goes extinct in the end.

Finally, we consider the stochastic model. the values of the parameters are same as those for Figure 2. And $\sigma_{1} = \sigma_{2} = \sigma_{3} = 0.1$. The following discretization equations in Milsteins type^[29] are used to iteratively calculate the approximate solutions of stochastic system (4.1) in Matlab programs:

where $\xi_{1, i}, \xi_{2, i},$ and $\xi_{3, i}$ are $N(0, 1)$-distributed independent Gaussian random variable, $\sigma_{1}, \sigma_{2}$ and $\sigma_{3}$ are intensities of white noise and time increment $\Delta t > 0$. The simulation results are shown in Figure 5.

6.   Conclusion

In this paper, firstly, phytoplankton-zooplankton model including TOM is proposed to investigate effects of TOM upon planktonic dynamics. By constructing Lyapunov functions and using LaSalle's Invariance Principle, global stability of the equilibria which is determined by threshold values $s_{0}$ and $s_{1}$ is established. The analytical results of $(2.2)$ indicate that TOM has significant effects on the lake ecosystem stability. Different constant input rates of TOM in the model can result in different results, which represents different biological meanings.

Plankton populations which reside in lake ecosystem are persistently influenced by the environmental fluctuations. In order to make more precise biological findings, we develop a stochastic model. We get that when the perturbation is small, the plankton populations reach a stationary distribution, which is ergodic. Then we simulate model (2.2) with many parameter values taken from the literature by varying the constant input rate of TOM along with the carrying capacity value. The simulations indicate that the magnitude of constant input rate of TOM and carrying capacity affect the persistence of phytoplankton populations greatly. We also simulate the stochastic model with the same parameter values. We can notice that mean value trajectory of (4.1) which is the average of ten trajectories oscillates randomly surrounding the solution of system (2.2), rather than reach the stable steady state (Figure 5). And the simulations also show that the trajectory of (4.1) is totally different from that of (2.2).

In summary, we set up exploratory models to consider the potential role of TOM that mediates interactions between trophic levels in a simple plankton food-chain. We have shown that, in principle, TOM plays an important role in influencing interactions between phytoplankton and zooplankton. The theoretical results of our deterministic model fit well with some experimental results. Further, when the perturbation is small, the stochastic model exists a stationary distribution which is ergodic. We emphasize that the qualitative behaviour of our models that we are interested in and this study are only meant as an initial exploratory attempt to consider the influence of TOM on the lake ecosystem.

Acknowledgments

This research is supported by the National Natural Sciences Foundation of China (Grant No. 11571326).

Conflict of interest

The authors declare there is no conflict of interest.