Mathematical model for the growth of Mycobacterium tuberculosis in the granuloma

1.   Introduction

Competitive relationships between different species are common among populations of organisms. Interspecific competition is a phenomenon in which two or more populations living together struggle for ecological resources. For example, crops and weeds in the same field are two populations competing with each other, and there is also competition between cattle and sheep in the same grassland. As a result, competition models are also important objects of study in mathematical biology. Mathematical models of interspecific competition can effectively explain and predict populations' changing rules and development trends in the process of mutual competition.

In the 1920s, Volterra proposed the famous Lotka-Volterra model based on the predator-prey relationship ^[1]. (Lotka had also proposed this model in his chemical reaction studies.) This model has resulted in a breakthrough in population model study and piqued the interest of a wide range of academics. After the predator-prey model, the Lotka-Volterra competition model was proposed and a series of studies on competition models were obtained. We refer the reader to ^{[2,3,4,5,6,7,8]} and the references therein. Specifically, a classic two-species competition model can be expressed as follows.

where ${{N}_{i}(t)}$ $(i = 1, 2)$ denotes the population size of the $i$th species at time $t$; ${{r}_{i}}$ is the intrinsic rate of increase; ${{b}_{ii}}$ denotes the intraspecific competition rate of the species. Since two populations compete with each other, ${{b}_{12}}$ denotes the effect of Population 2 on Population 1, and ${{b}_{21}}$ denotes the effect of Population 1 on Population 2. In general, ${{b}_{12}}\ne {{b}_{21}}$. All parameters in Eq (1.1) are positive. Regarding the equilibrium points of Eq (1.1), we have the following conclusions ^[9].

1) Equilibrium ${{E}_{0}} = (0, 0)$ always exists.

2) Equilibrium ${{E}_{1}} = (\frac{{{r}_{1}}}{{{b}_{11}}}, 0)$ always exists.

3) Equilibrium ${{E}_{2}} = (0, \frac{{{r}_{2}}}{{{b}_{22}}})$ always exists.

4) Define $\alpha = {{b}_{22}}{{r}_{1}}-{{b}_{12}}{{r}_{2}}$, $\beta = -{{b}_{21}}{{r}_{1}}+{{b}_{11}}{{r}_{2}}$ and $\gamma = {{b}_{11}}{{b}_{22}}-{{b}_{12}}{{b}_{21}}$. If the condition $\alpha > 0$, $\beta > 0$ or $\alpha < 0$, $\beta < 0$ holds, equilibrium ${{E}^{*}} = (N_{1}^{*}, N_{2}^{*}) = (\frac{\alpha }{\gamma }, \frac{\beta }{\gamma })$ always exists.

The discussion of Eq (1.1) does not consider the effects of environmental disturbances, while in real life, the dynamic behaviors of biological populations are inevitably affected by environmental disturbances such as air humidity and weather changes. Several scholars' studies have proved this view. Mode and Jacobson ^[10] found that the probability of population extinction is sensitive to environmental perturbations. DuBowy ^[11] found that seasonal shifts affect population growth and community composition. And Mao et al. ^[12] showed that environmental noise can inhibit population size surges. Therefore, the introduction of stochastic factors is essential.

There are mainly two ways to characterize environmental disturbances. One method is to directly introduce environmental noises to a deterministic model, and it has been studied and used by many scholars ^{[13,14,15,16,17,18]}. Another approach is to simulate random environmental perturbations through the use of a stochastic process. Allen's study ^[19] illustrated that a stochastic Ornstein-Uhlenbeck process (also known as the mean-reverting process) possesses several advantages over linear functions of white noise in terms of being able to modify parameters for environmental variability, which indicated that the second way is a feasible and biologically meaningful approach. The study by Zhang and Yuan ^[20] showed that the Ornstein-Uhlenbeck process is an effective and reasonable way to introduce environmental noise into the continuous culture model of microorganisms. It was also found that the reversion speed and volatility intensity have essential effects on the extinction and persistence of microorganisms. Song and Zhang ^[21] proposed a new stochastic SVEIS model with the Ornstein-Uhlenbeck process and studied this model's stationary distribution and extinction. Yang et al. ^[22] introduced the Ornstein-Uhlenbeck process into a food chain system to simulate stochastic perturbation, and some valuable theoretical results were obtained. For the related study on the Ornstein-Uhlenbeck process, readers can refer to ^{[19,20,21,22,23,24,25,26,27,28]}. Previous studies have shown that introducing the Ornstein-Uhlenbeck process into the model is worth investigating in depth. However, as far as we know, no one has studied a competition model with the Ornstein-Uhlenbeck process.

In fact, the growth rate and death rate of the population are easily disturbed by environmental changes. Therefore, we consider a comparison of two methods to incorporate environmental variability into the intrinsic rates of increase. There are some differences between the Ornstein-Uhlenbeck noise and conventional Gaussian white noise, and these differences made us prefer to use the Ornstein-Uhlenbeck process to model this environmental variability.

1) The first approach is to introduce Gaussian white noise, specifically by making the following transformation to ${{r}_{1}}$ and ${{r}_{2}}$:

where ${{B}_{r1}(t)}$ and ${{B}_{r2}(t)}$ are mutually independent Brownian motions. By directly integrating Eq (1.2), the average per intrinsic rates of increase over an interval $[0, T]$ is equal to

That is, the variances of the average intrinsic rates of increase ${{\bar{r}}_{1}}$ and ${{\bar{r}}_{2}}$ over an interval of length $T$ both tend to infinity as $T\to 0$.

Specifically operating the second method, i.e., we assume that ${{r}_{1}}$ and ${{r}_{2}}$ have a transformation of the following form:

${{m}_{1}(t)}$ and ${{m}_{2}(t)}$ are the Ornstein-Uhlenbeck processes which satisfy

where ${{\theta }_{i}}$ is the speed of reversion and ${{\xi }_{i}}$ is the intensity of volatility of the process ${{m}_{i}}(t)$ (${{\theta }_{i}} > 0, {{\xi }_{i}} > 0, i = 1, 2$). ${{B}_{1}}(t)$ and ${{B}_{2}}(t)$ are mutually independent Brownian motions. After calculation, Eq (1.4) can be solved exactly to yield

where ${{m}_{i0}}$ is the initial value of the process ${{m}_{i}}(t)$. From Eq (1.5), we obtain that ${{m}_{i}}(t)\sim N({{m}_{i0}}{{e}^{-{{\theta }_{i}}t}}$, $\frac{{{\theta }_{i}}^{2}}{2{{\xi }_{i}}}(1-{{e}^{-2{{\theta }_{i}}t}}))$. It is not difficult to find that the process ${{m}_{i}}(t)$ follows the distribution $N(0, \frac{{{\theta }_{i}}^{2}}{2{{\xi }_{i}}})$ as $t$ tends to infinity. By directly integrating Eq (1.5), the average per intrinsic rates of increase over an interval $[0, T]$ is equal to

For small values of $T$, we have

It is not difficult to find that, unlike the results for Gaussian white noise, the variance goes to 0 rather than $\infty$ as $T\to 0$.

2) For a small time interval $\Delta t$, an Ornstein-Uhlenbeck process $X(t)$ has a correlation coefficient $\rho (X(t), X(t+\Delta t)) = 1-o(\Delta t)$ while $\rho (X(t), X(t+\Delta t)) = 0$ for Gaussian white noise. So white noise processes are often used to model random perturbations with very small correlation periods. However, the environment of many biological systems can be considered as continuously varying due to a large number of interacting variables. Therefore, the Ornstein-Uhlenbeck process is more suitable for modeling the parameters in this paper.

Motivated by the above discussion, we consider using the Ornstein-Uhlenbeck process to simulate the stochastic process. The interspecific competition model with the Ornstein-Uhlenbeck process studied in this paper is obtained by combining Eqs (1.1), (1.3) and (1.4):

Throughout this paper, we let $(\Omega, F, {{\left\{ {{F}_{t}} \right\}}_{t\ge 0}}, P)$ be a complete probability space with a filtration ${{\left\{ {{F}_{t}} \right\}}_{t\ge 0}}$ satisfying the usual conditions (i.e., it is right continuous and increasing while ${{F}_{0}}$ contains all $P$-null sets). And we denote $R_{+}^{n} = \left\{ x\in {{R}^{n}}:{{x}_{i}} > 0, 1\le i\le n \right\}$, ${{R}_{+}} = [0, +\infty)$.

This paper is organized as follows. In the next section, we investigate the existence and uniqueness of the global solution to the model (1.6). In Section 3, we discuss the cases of population extinction. In Section 4, we obtain sufficient conditions for the existence of the stationary distribution; the expression for the mean and covariance of the density function of the linearized system corresponding to the stochastic model (1.6) around the original point are obtained in Section 5. In Section 6, the previous theoretical results are verified by numerical simulations. Finally, a brief conclusion is given.

2.   Existence and uniqueness of the global positive solution

Since ${{N}_{1}}(t)$ and ${{N}_{2}}(t)$ denote the number of individuals, they should be non-negative from the viewpoint of biology. To study the long-term dynamical behaviors of the stochastic model (1.6) proposed in this paper, we first need to consider whether the solution is global and non-negative, which is a fundamental condition for the follow-up study. So in this section, we obtain the following theorem, which guarantees the existence and uniqueness of the global positive solution for the model (1.6).

Theorem 1. For any given initial value $({{N}_{1}}(0), {{N}_{2}}(0), {{m}_{1}}(0), {{m}_{2}}(0))\in R_{+}^{2}\times {{R}^{2}}$, the model (1.6) has a unique solution $({{N}_{1}}(t), {{N}_{2}}(t), {{m}_{1}}(t), {{m}_{2}}(t))\in R_{+}^{2}\times {{R}^{2}}$ for all $t\ge 0$ with a probability of one.

Proof. Apparently, coefficients of the model (1.6) satisfy the local Lipschitz condition; so, for any given initial value $({{N}_{1}}(0), {{N}_{2}}(0), {{m}_{1}}(0), {{m}_{2}}(0))\in R_{+}^{2}\times {{R}^{2}}$, there exists a unique solution $({{N}_{1}}(t), {{N}_{2}}(t), {{m}_{1}}(t), {{m}_{2}}(t))\in R_{+}^{2}\times {{R}^{2}}$, $t\in [0, {{\rho }_{e}})$, where ${{\rho }_{e}}$ is the explosion time. To prove the global nature of the solution, we only need to prove that ${{\rho }_{e}} = \infty \text{ }a.s.$ Let ${{l}_{0}} > 0$ be sufficiently large such that ${{N}_{1}}(0), {{N}_{2}}(0), {{e}^{{{m}_{1}}(0)}}$ and ${{e}^{{{m}_{2}}(0)}}$ all lie within the interval $[\frac{1}{{{l}_{0}}}, {{l}_{0}}]$. For each integer $l > {{l}_{0}}$, define the stopping time ^[29] as follows:

where $\inf \phi = \infty$. If we can show that ${{\tau }_{\infty }} = \infty \text{ }a.s.$, then ${{\rho }_{e}} = \infty \text{ }a.s.$, and, further, Theorem 1 is proved.

Now, considering the contradiction, we assume that there exists a pair of constants $T > 0$ and $\varepsilon \in (0, 1)$ such that $P({{\tau }_{\infty }}\le T) > \varepsilon$. Hence there is an integer ${{l}_{1}}\ge {{l}_{0}}$ such that $P({{\tau }_{l}}\le T)\ge \varepsilon$, $\forall l\ge {{l}_{1}}$. Define ${{C}^{2}}$-function $V:R_{+}^{2}\times {{R}^{2}}\to {{R}_{+}}$,

Applying Itô's formula to V, we get

where

${{k}_{1}}$ is a positive constant that is independent of the initial value. Integrating both sides of Eq (2.1) from 0 to ${{\tau }_{l}}\wedge T$ and then taking expectations, we get

Let $\tilde{\Omega } = \left\{ \varpi \in \Omega :{{\tau }_{l}}\le T \right\}$, where $\varpi$ represents a sample point; then, we have $P(\tilde{\Omega })\ge \varepsilon$, $\varepsilon \in (0, 1)$; thus,

where ${{I}_{{\tilde{\Omega }}}}$ is the indicator function of $\tilde{\Omega }$. Let $l\to \infty$; we obtain

which is a contradiction. This completes the proof of Theorem 1.

Theorem 2. For any given initial value $({{N}_{1}}(0), {{N}_{2}}(0), {{m}_{1}}(0), {{m}_{2}}(0))\in R_{+}^{2}\times {{R}^{2}}$, the solution $({{N}_{1}}(t), {{N}_{2}}(t)$, ${{m}_{1}}(t), {{m}_{2}}(t))$ of the model (1.6) has the property that for any $p\ge 1$, there exists a constant $k(p)$ such that

Moreover,

Proof. Define a Lyapunov function as

where $p\ge 1$. Applying Itô's formula to V, we get

where

For any constant $\delta$ which satisfies $0 < \frac{\delta }{4p} < \frac{{{\theta }_{i}}}{2}, i = 1, 2$, we have

where

Thus,

which implies

There exists a constant $h(p) > 0$ such that $E[{{({{N}_{1}(t)}+{{N}_{2}(t)})}^{p}}]\le h(p), t\ge 0\text{ }a.s.$ For convenience, define $N(t) = {{({{N}_{1}}(t)+{{N}_{2}}(t))}^{p}} = {{x}^{p}(t)}$. Applying Itô's formula to $N(t)$ yields

Let $\theta > 0$ be sufficiently small and satisfy $n\theta \le t\le (n+1)\theta, n = 1, 2, \ldots$. It follows that

where

Choose $\theta$ sufficiently small such that $I < h(p)$; therefore,

Let $\varepsilon$ be an arbitrary positive constant; then, based on Chebyshev's inequality ^[30], it follows that

There exists an integer-valued random variable ${{n}_{0}}(\omega)$ such that for almost all $\omega \in \Omega$, when $n\ge {{n}_{0}}$, we have

If $n\ge {{n}_{0}}$ and $n\theta \le t\le (n+1)\theta$, we get

Let $\varepsilon \to 0$; we get

then

Thus,

and it follows that

3.   Extinction of the model (1.6)

In this section, the extinction of populations is discussed. First, we give the following lemma needed for the subsequent proof.

Lemma 1. ^[31] If $f(t)$ is integrable, we define ${{\langle f \rangle }_{t}} = \frac{1}{t}\int_{0}^{t}{f(s)\text{d}s}$. Assume that $z(t)\in C(\Omega \times [0, \infty), {{R}_{+}})$.

1) If for all $t > T$, there are two positive constants $T$ and ${{\delta }_{0}}$ satisfying

where ${{\alpha }_{i}}$ and $\delta$ are constants, we get

2) If for all $t > T$, there are three positive constants $T$, $\delta$ and ${{\delta }_{0}}$ satisfying

where ${{\alpha }_{i}}$ denotes constants, we get $\underset{t\to \infty }{\mathop{\lim \inf }}\, {{\langle z \rangle }_{t}}\ge \frac{\delta }{{{\delta }_{0}}}\text{ }a.s.$

Theorem 3. For any given initial value $({{N}_{1}}(0), {{N}_{2}}(0), {{m}_{1}}(0), {{m}_{2}}(0))\in R_{+}^{2}\times {{R}^{2}}$, the model (1.6) has the following properties:

1) if $\alpha > 0$ and $\beta < 0$, then

2) if $\alpha < 0$ and $\beta > 0$, then

Proof. 1) The discussion is divided into the following two cases.

Case 3.1 If $\alpha > 0$, $\beta < 0$ and $\gamma > 0$, by Itô's formula, we get

According to Theorem 2, let ${{\varepsilon }_{1}}$ satisfy $\frac{-\beta }{{{b}_{11}}} = \frac{{{b}_{21}}{{r}_{1}}-{{b}_{11}}{{r}_{2}}}{{{b}_{11}}} > {{\varepsilon }_{1}} > 0$ and there exist ${{T}_{1}}$ which is sufficiently large such that $\frac{{{b}_{21}}}{{{b}_{11}}t}\ln \frac{{{N}_{1}}(t)}{{{N}_{1}}(0)}+\frac{1}{t}\ln {{N}_{2}}(0) < {{\varepsilon }_{1}}$ holds for $t > {{T}_{1}}$. It is easy to show that

Taking the superior limit on both sides, we have

Therefore, we have that $\underset{t\to \infty }{\mathop{\lim }}\, {{N}_{2}}(t) = 0\text{ }a.s.$ Considering the equation

according to Lemma 1, we have

Case 3.2 Similar to Case 3.1, if $\alpha > 0$, $\beta < 0$ and $\gamma > 0$, by Itô's formula, we get

According to Theorem 2, let ${{\varepsilon }_{2}}$ satisfy $\frac{\alpha }{{{b}_{22}}} = \frac{{{b}_{22}}{{r}_{1}}-{{b}_{12}}{{r}_{2}}}{{{b}_{22}}} > {{\varepsilon }_{2}} > 0$ and there exist ${{T}_{2}}$ which is sufficiently large such that $\frac{1}{t}\ln \frac{{{N}_{1}}(t)}{{{N}_{1}}(0)}+\frac{{{b}_{12}}}{{{b}_{22}}t}\ln {{N}_{2}}(0) < {{\varepsilon }_{2}}$ holds for $t > {{T}_{2}}$. It is easy to show that

Taking the superior limit on both sides, we have

Therefore, we have $\mathop {\lim }\limits_{t \to \infty } {\mkern 1mu} {{N}_{2}}(t) = 0\text{ }a.s.$ Furthermore, we have

2) The discussion is divided into the following two cases.

Case 3.3 If $\alpha < 0$, $\beta > 0$ and $\gamma > 0$, by Itô's formula, we get

According to Theorem 2, let ${{\varepsilon }_{3}}$ satisfy $\frac{-\alpha }{{{b}_{12}}} = \frac{-{{b}_{22}}{{r}_{1}}+{{b}_{12}}{{r}_{2}}}{{{b}_{12}}} > {{\varepsilon }_{3}} > 0$ and there exist ${{T}_{3}}$ which is sufficiently large such that $\frac{{{b}_{22}}}{{{b}_{12}}t}\ln {{N}_{1}}(0)+\frac{1}{t}\ln \frac{{{N}_{2}}(t)}{{{N}_{2}}(0)} < {{\varepsilon }_{3}}$ holds for $t > {{T}_{3}}$. It is easy to show that

Taking the superior limit on both sides, we have

Therefore, we have that $\mathop {\lim }\limits_{t \to \infty } {\mkern 1mu} {{N}_{1}}(t) = 0\text{ }a.s.$ Considering the equation

according to Lemma 1, we have

Case 3.4 Similar to Case 3.3, if $\alpha < 0$, $\beta > 0$ and $\gamma < 0$, by Itô's formula, we get

According to Theorem 2, let ${{\varepsilon }_{4}}$ satisfy $\frac{\beta }{{{b}_{21}}} = \frac{-{{b}_{21}}{{r}_{1}}+{{b}_{11}}{{r}_{2}}}{{{b}_{21}}} > {{\varepsilon }_{4}} > 0$ and there exist ${{T}_{4}}$ which is sufficiently large such that $\frac{1}{t}\ln {{N}_{1}}(0)+\frac{{{b}_{11}}}{{{b}_{21}}t}\ln \frac{{{N}_{2}}(t)}{{{N}_{2}}(0)} < {{\varepsilon }_{4}}$ holds for $t > {{T}_{4}}$. It is easy to show that

Taking the superior limit on both sides, we have

Therefore, we have that $\mathop {\lim }\limits_{t \to \infty } {\mkern 1mu} {{N}_{1}}(t) = 0\text{ }a.s.$ Furthermore, we have

4.   Stationary distribution of the model (1.6)

In this section, a sufficient condition for the existence of the stationary distribution in the model (1.6) is obtained. Before giving the theorem, we present the following lemma needed for the subsequent proof.

Lemma 2. ^[32] $X(t)$ is a homogeneous Markov process which is expressed as

Assume that $b(X), {{\sigma }_{1}}(X), \cdots, {{\sigma }_{k}}(X)$$(t\ge {{t}_{0}}, x\in {{R}^{d}})$ are continuous.

1) There is a constant $B$ that satisfies

2) A non-negative $U(x)\in {{C}^{2}}$ in ${{R}^{d}}$ that satisfies $LU(x)\le -1$ outside of some compact set exists.

If Conditions (1) and (2) hold, $X(t)$ is a stationary Markov process.

Theorem 4. For any given initial value $({{N}_{1}}(0), {{N}_{2}}(0), {{m}_{1}}(0), {{m}_{2}}(0))\in R_{+}^{2}\times {{R}^{2}}$, suppose that the condition $\alpha > 0, \beta > 0$ holds. If

where $w = \frac{{{\xi }_{1}}}{2{{\theta }_{1}}}+\frac{{{\xi }_{2}}}{2{{\theta }_{2}}}$, then the model (1.6) has a stationary distribution.

Proof. For the model (1.6), Condition (1) of Lemma 2 holds. Thus, we only need to verify Condition (2). Consider the following function:

where $N_{1}^{*} = \frac{{{b}_{22}}{{r}_{1}}-{{b}_{12}}{{r}_{2}}}{{{b}_{11}}{{b}_{22}}-{{b}_{12}}{{b}_{21}}} > 0$ and $N_{2}^{*} = \frac{-{{b}_{21}}{{r}_{1}}+{{b}_{11}}{{r}_{2}}}{{{b}_{11}}{{b}_{22}}-{{b}_{12}}{{b}_{21}}} > 0$. For the sake of convenience, define

Using Itô's formula, we get

Therefore,

where $w = \frac{{{\xi }_{1}}}{2{{\theta }_{1}}}+\frac{{{\xi }_{2}}}{2{{\theta }_{2}}}$. If

holds, then the ellipsoid $-({{b}_{11}}-\frac{{{b}_{12}}+{{b}_{21}}}{2}-\frac{{{\xi }_{1}}}{2}){{({{N}_{1}}-N_{1}^{*})}^{2}}-({{b}_{22}}-\frac{{{b}_{12}}+{{b}_{21}}}{2}-\frac{{{\xi }_{2}}}{2}){{({{N}_{2}}-N_{2}^{*})}^{2}}-\frac{1}{2{{\xi }_{1}}}m_{1}^{2}-\frac{1}{2{{\xi }_{2}}}m_{2}^{2}+w$ lies entirely in $R_{_{+}}^{2}\times {{R}^{2}}$. We can take $U$ to be a neighborhood of the ellipsoid with $\bar{U}\subseteq R_{+}^{2}\times {{R}^{2}}$ where $\bar{U}$ represents the compact closure of $U$. Model (1.6) is known to satisfy Condition (1) of Lemma 2, and we have just proved that Condition (2) is true. Hence, the model (1.6) has a stationary distribution according to Lemma 2.

5.   Density function of the model (1.6)

In this section, the formulas for the mean and the covariance of the probability density function of the corresponding linearized system near the equilibrium point are obtained. For the model (1.6), the equilibrium is ${{E}^{*}} = (N_{1}^{*}, N_{2}^{*}, 0, 0)$. Let ${{u}_{i}} = {{N}_{i}}-N_{i}^{*}, i = 1, 2$; the corresponding linearized model of the model (1.6) is

where ${{a}_{11}} = {{b}_{11}}N_{1}^{*} > 0$, ${{a}_{12}} = {{b}_{12}}N_{1}^{*} > 0$, ${{a}_{21}} = {{b}_{21}}N_{2}^{*} > 0$ and ${{a}_{22}} = {{b}_{22}}N_{2}^{*} > 0$.

Theorem 5. For any given initial value $({{N}_{1}}(0), {{N}_{2}}(0), {{m}_{1}}(0), {{m}_{2}}(0))\in R_{+}^{2}\times {{R}^{2}}$, if the conditions of Theorem 4 hold, then the model (5.1) has a stationary distribution around the original point. The mean vector is $(0, 0, 0, 0)$, and the covariance matrix has the following form:

where

${{p}_{1}} = {{a}_{11}}+{{a}_{22}}$, ${{p}_{2}} = {{a}_{11}}{{a}_{22}}-{{a}_{21}}{{a}_{12}}$, ${{\alpha }_{1}} = {{a}_{21}}{{\xi }_{1}}N_{1}^{*}$, ${{\alpha }_{2}} = {{a}_{12}}{{\xi }_{2}}N_{2}^{*}$, ${{\eta }_{1}}\text{ = }2\theta _{1}^{2}{{p}_{1}}+2{{\theta }_{1}}p_{1}^{2}+2{{p}_{1}}{{p}_{2}}$, ${{\eta }_{2}}\text{ = }2\theta _{2}^{2}{{p}_{1}}+2{{\theta }_{1}}p_{1}^{2}+2{{p}_{1}}{{p}_{2}}$,

Proof. Let $X = {{({{u}_{1}}, {{u}_{2}}, {{m}_{1}}, {{m}_{2}})}^{T}}$, $B(t) = {{(0, 0, {{B}_{1}}(t), {{B}_{2}}(t))}^{T}}$ and

then the model (5.1) is transformed into $\text{d}X = AX(t)\text{d}t+H\text{d}B(t)$. Then we get

Since ${{B}_{1}}(t)$ and ${{B}_{2}}(t)$ are mutually independent, Equation (5.2) can be equivalently transformed into the following algebraic sub-equations:

where $H_{1}^{2} = diag(0, 0, \xi _{1}^{2}, 0)$, $H_{2}^{2} = diag(0, 0, 0, \xi _{2}^{2})$ and $\Sigma = {{\Sigma }_{1}}+{{\Sigma }_{2}}$. Let $B = \left(\begin{matrix} -{{a}_{11}} & -{{a}_{12}} \\ -{{a}_{21}} & -{{a}_{22}} \\ \end{matrix} \right)$; the characteristic polynomial of matrix $B$ is

where ${{p}_{1}} = {{a}_{11}}+{{a}_{22}} > 0$, ${{p}_{2}} = {{a}_{11}}{{a}_{22}}-{{a}_{12}}{{a}_{21}} > 0$ and ${{p}_{1}}{{p}_{2}} = ({{a}_{11}}+{{a}_{22}})({{a}_{11}}{{a}_{22}}-{{a}_{12}}{{a}_{21}}) > 0$. According to the Routh-Hurwitz criterion ^[33], matrix $B$ has all negative real-part eigenvalues. Next, we prove Theorem 5 in two cases.

Case 5.1 Consider the equation

Let ${{A}_{1}} = {{I}_{1}}AI_{1}^{-1}$, where ${{I}_{1}} = \left(\begin{matrix} 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right)$; we have

Next, let ${{A}_{2}} = {{I}_{2}}{{A}_{1}}I_{2}^{-1}$, where ${{I}_{2}} = \left(\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & -{{a}_{21}} & -{{a}_{22}} & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right)$; we therefore have

Let ${{A}_{3}} = {{I}_{3}}{{A}_{2}}I_{3}^{-1}$, where ${{I}_{3}} = \left(\begin{matrix} -{{a}_{21}}N_{1}^{*} & -{{p}_{1}} & -{{p}_{2}} & -{{a}_{22}}N_{2}^{*} \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right)$; thus,

Equation (5.3) is transformed into

Denote ${{\Sigma }_{01}}\text{ = }\frac{1}{\alpha _{1}^{2}}({{I}_{3}}{{I}_{2}}{{I}_{1}}){{\Sigma }_{1}}{{({{I}_{3}}{{I}_{2}}{{I}_{1}})}^{T}}$, where ${{\alpha }_{1}} = {{a}_{21}}{{\xi }_{1}}N_{1}^{*}$; we can get

where ${{\eta }_{1}}\text{ = }2\theta _{1}^{2}{{p}_{1}}+2{{\theta }_{1}}p_{1}^{2}+2{{p}_{1}}{{p}_{2}}$.

${{\Sigma }_{01}}$ is a positive semi-definite matrix, and ${{\Sigma }_{1}}\text{ = }\alpha _{1}^{2}{{({{I}_{3}}{{I}_{2}}{{I}_{1}})}^{-1}}{{\Sigma }_{01}}{{[{{({{I}_{3}}{{I}_{2}}{{I}_{1}})}^{-1}}]}^{T}}$ is also a positive semi-definite matrix.

Case 5.2 Consider the equation

Let ${{B}_{1}} = {{J}_{1}}AJ_{1}^{-1}$, where ${{J}_{1}} = \left(\begin{matrix} 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \end{matrix} \right)$; we have

Next, let ${{B}_{2}} = {{J}_{2}}{{B}_{1}}J_{2}^{-1}$, where ${{J}_{2}} = \left(\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & -{{a}_{12}} & -{{a}_{11}} & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right)$; we therefore have

Let ${{B}_{3}} = {{J}_{3}}{{B}_{2}}J_{3}^{-1}$, where ${{J}_{3}} = \left(\begin{matrix} -{{a}_{12}}N_{2}^{*} & -{{p}_{1}} & -{{p}_{2}} & -{{a}_{11}}N_{1}^{*} \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right)$; thus,

Equation (5.4) is transformed into

Denote ${{\Sigma }_{02}}\text{ = }\frac{1}{\alpha _{2}^{2}}({{J}_{3}}{{J}_{2}}{{J}_{1}}){{\Sigma }_{2}}{{({{J}_{3}}{{J}_{2}}{{J}_{1}})}^{T}}$, where ${{\alpha }_{2}} = {{a}_{12}}{{\xi }_{2}}N_{2}^{*}$; we can get

where ${{\eta }_{2}}\text{ = }2\theta _{2}^{2}{{p}_{1}}+2{{\theta }_{1}}p_{1}^{2}+2{{p}_{1}}{{p}_{2}}$.

${{\Sigma }_{02}}$ is a positive semi-definite matrix, and ${{\Sigma }_{2}}\text{ = }\alpha _{2}^{2}{{({{J}_{3}}{{J}_{2}}{{J}_{1}})}^{-1}}{{\Sigma }_{02}}{{[{{({{J}_{3}}{{J}_{2}}{{J}_{1}})}^{-1}}]}^{T}}$ is also a positive semi-definite matrix. The proof is complete.

6.   Numerical simulations

This section gives examples of numerical simulations to verify the previous results. Specifically, we will verify the following results:

1) The effect of the change of ${{\xi }_{i}}$ in the Ornstein-Uhlenbeck process.

2) When the conditions of Theorems 3 or 4 are satisfied, observe whether the images resulting from numerical simulation agree with the theory.

3) Compare and observe the images of the deterministic and stochastic models.

The corresponding discrete model of the model (1.6) is

where $\Delta t > 0$ is the time increment and ${{\nu }_{k}}$ and ${{l}_{k}}$ are the mutually independent Gaussian random variables with distribution $N(0, 1)$ for $k = 1, 2, \ldots, n$. $N_{1}^{k}, N_{2}^{k},$ $m_{1}^{k}$ and $m_{2}^{k}$ denote the corresponding values of the $k$th iteration of the model (6.1).

And the corresponding discrete model of the deterministic model (1.1) is

where $\Delta t > 0$ is the time increment and $N_{1}^{k}$ and $N_{2}^{k}$ denote the corresponding values of the $k$th iteration of the model (6.2).

1) We will verify the correctness of Theorem 4 and consider the effect of the change of ${{\xi }_{i}}$ in the Ornstein-Uhlenbeck process.

Example 1. Let ${{r}_{1}} = 0.4, {{r}_{2}} = 0.5, {{b}_{11}} = 0.5, {{b}_{12}} = 0.4, {{b}_{21}} = 0.45, {{b}_{22}} = 0.6, {{\theta }_{1}} = {{\theta }_{2}} = 0.5, {{N}_{1}}(0) = 0.5, {{N}_{2}}(0) = 0.5, {{M}_{1}}(0) = 0.001$ and ${{M}_{2}}(0) = 0.001$. If ${{\xi }_{1}} = {{\xi }_{2}} = 0.001$, then $\alpha = 0.04 > 0,$ $\beta = 0.07 > 0$ and $w = 0.002 < \min \{0.0083, 0.0594\}$. The corresponding image can be seen in Figure 1.

Example 2. Similarly, let ${{r}_{1}} = 0.4, {{r}_{2}} = 0.5, {{b}_{11}} = 0.5, {{b}_{12}} = 0.4, {{b}_{21}} = 0.45, {{b}_{22}} = 0.6, {{\theta }_{1}} = {{\theta }_{2}} = 0.5, {{N}_{1}}(0) = 0.5, {{N}_{2}}(0) = 0.5, {{M}_{1}}(0) = 0.001$ and ${{M}_{2}}(0) = 0.001$. If ${{\xi }_{1}} = 0.005$ and ${{\xi }_{2}} = 0.001$, then $\alpha = 0.04 > 0,$ $\beta = 0.07 > 0$ and $w = 0.006 < \min\left\{0.0081, 0.0594 \right\}$. The corresponding image can be seen in Figure 2.

Example 3. Similarly, let ${{r}_{1}} = 0.4, {{r}_{2}} = 0.5, {{b}_{11}} = 0.5, {{b}_{12}} = 0.4, {{b}_{21}} = 0.45, {{b}_{22}} = 0.6, {{\theta }_{1}} = {{\theta }_{2}} = 0.5, {{N}_{1}}(0) = 0.5, {{N}_{2}}(0) = 0.5, {{M}_{1}}(0) = 0.001$ and ${{M}_{2}}(0) = 0.001$. If ${{\xi }_{1}} = {{\xi }_{2}} = 0.005$, then $\alpha = 0.04 > 0,$ $\beta = 0.07 > 0$ and $w = 0.01 > \min \{0.0081, 0.0587\}$. The corresponding image can be seen in Figure 3.

Figures 1 and 2 show that if the conditions of Theorem 4 hold, the two species of the model (1.6) can coexist in the long term. By comparing Figures 1–3, we find that when other parameters and initial values remain the same, the intensity of the fluctuations in the plot increases if ${{\xi }_{i}}$, i.e., is the intensity of the volatility of the process increases.

2) There are two examples of extinction.

Example 4. Let ${{r}_{1}} = 0.4$, ${{r}_{2}} = 0.5$, ${{b}_{11}} = 0.5$, ${{b}_{12}} = 0.4$, ${{b}_{21}} = 0.7$, ${{b}_{22}} = 0.6$, ${{\theta }_{1}} = {{\theta }_{2}} = 0.5$, ${{\xi }_{1}} = {{\xi }_{2}} = 0.01$, ${{N}_{1}}(0) = 0.5$, ${{N}_{2}}(0) = 0.5$, ${{M}_{1}}(0) = 0.001$ and ${{M}_{2}}(0) = 0.001$; then, $\alpha = 0.04 > 0$ and $\beta = -0.03 < 0$.

By Theorem 3, Population 2 will go extinct eventually and

The corresponding image can be seen in Figure 4.

Example 5. Let ${{r}_{1}} = 0.4$, ${{r}_{2}} = 0.5$, ${{b}_{11}} = 0.5$, ${{b}_{12}} = 0.5$, ${{b}_{21}} = 0.61$, ${{b}_{22}} = 0.6$, ${{\theta }_{1}} = {{\theta }_{2}} = 0.5$, ${{\xi }_{1}} = {{\xi }_{2}} = 0.01$, ${{N}_{1}}(0) = 0.5$, ${{N}_{2}}(0) = 0.5$, ${{M}_{1}}(0) = 0.001$ and ${{M}_{2}}(0) = 0.001$; then, $\alpha = -0.01 < 0$ and $\beta = 0.01 > 0$.

By Theorem 3, Population 1 will go extinct eventually and

The corresponding image can be seen in Figure 5.

From Figures 4 and 5, it is not difficult to find that the result is in line with Theorem 3.

3) We consider the comparison of the deterministic model (1.1) and the stochastic model (1.6).

Example 6. Consider the discrete models (6.1) and (6.2), take the parameters as ${{r}_{1}} = 0.4$, ${{r}_{2}} = 0.5$, ${{b}_{11}} = 0.5$, ${{b}_{12}} = 0.4$, ${{b}_{21}} = 0.45$, ${{b}_{22}} = 0.6$, ${{\theta }_{1}} = {{\theta }_{2}} = 0.5$ and ${{\xi }_{1}} = {{\xi }_{2}} = 0.0006$ and set the initial values as ${{N}_{1}}(0) = 0.5$, ${{N}_{2}}(0) = 0.5$, ${{M}_{1}}(0) = 0.001$ and ${{M}_{2}}(0) = 0.001.$ The corresponding image can be seen in Figure 6.

It is not difficult to find that when the intensity of the volatility of the Ornstein-Uhlenbeck process is relatively small, the trends of the stochastic and deterministic model populations are similar.

7.   Conclusions

Competition models are an important medium for us to study the extinction and survival of competing species in ecosystems, and they have theoretical support for the study of ecosystem balance. This paper studies the dynamical behaviors of two population competition models with the Ornstein-Uhlenbeck process. The existence and uniqueness of the system solution have been verified, the cases of population extinction have been discussed and sufficient conditions for the stationary distribution of the system have been obtained. Specifically, if $\alpha > 0$ and $\beta < 0$ or $\alpha < 0$ and $\beta > 0$, then one population goes extinct while the other survives. Otherwise, if $\alpha > 0$, $\beta > 0$ and the parameters ${{\theta }_{i}}$ and ${{\xi }_{i}}$ of the Ornstein-Uhlenbeck process meet certain conditions, then the model (1.6) exists as a stable distribution, which means that the two populations coexist in this case.

In addition to proving the existence and uniqueness of global positive solutions and obtaining sufficient conditions for population extinction and stationary distribution, the innovation of this paper is that we further obtained the results of the mathematical analysis of the density function of the stochastic model (1.6), which is challenging work. Finally, the numerical simulations with examples have been provided to verify the theoretical results of extinction and stationary distribution. It is also proved that the stochastic and the corresponding deterministic models have similar properties when the intensity of the volatility of the Ornstein-Uhlenbeck process is relatively small.

Previously, although using the Ornstein-Uhlenbeck process to model stochastic processes is an effective way to introduce environmental noise into stochastic models, there are fewer related studies. Moreover, no one has studied the dynamical behaviors of stochastic competition models with Ornstein-Uhlenbeck processes, so this paper has some reference significance for studying stochastic competition models.

Acknowledgments

Thanks to the Reviewers for their helpful comments and suggestions. This work was supported by the National Natural Science Foundation of China Tianyuan Mathematical Foundation (No. 12126312, No. 12126328). The authors gratefully acknowledge the Natural Science Foundation of Heilongjiang Province (No. LH2022E023) for the support in publishing this paper.

Conflict of interest

The authors declare that there is no conflict of interest.