Dynamics of two phytoplankton populations under predation

1.   Introduction

Epidemic refers to the disease that can spread in a certain scale within a relatively short period of time. It caused a major hazard to human health. For example, the H1N1 influenza virus, which erupted in Mexico in 2009 and soon spread around the globe, is a dangerous virus and can develop to pneumonia. Some infectives may even appear respiratory failure, multiple organ damage, and finally be kicked to death. The H1N1 has killed about 16 thousand people throughout the world then and over 1.3 million people had been infected. Studying the transmission mechanism of epidemics can help reduce the spread of the diseases. Many researchers use the mathematical models to describe the transmission of epidemics ^[1], and various methods are derived to investigate the dynamics of the diseases ^[2,3].

Researchers have proposed mathematical models of infectious diseases and studied their dynamical behaviours extensively ^[4,5]. The dynamical behaviours of ordinary differential equations have been widely studied, and abundant dynamical behaviours have been found, such as periodic oscillations, bifurcations, stable limit cycles and time-delay effects ^[6,7,8,9,10]. To investigate the dynamics of cholera which spread in the European Mediterranean regions, Capasso and Paveri-fontana proposed an epidemic model in 1979 as follows ^[6]:

where $ u $ and $ v $ represent the concentration of the infectious agent and the infected population, respectively. $ \frac{1}{\alpha_1} $ denotes the average lifetime of the infectious agent, $ \frac{1}{\alpha_2} $ denotes the average infectious period of infected people, thus $ \alpha_1 $, $ \alpha_2 > 0 $. $ a $ denotes the multiplicative factor of the virus caused by the increasing number of infected people. The function $ g(u) $ represents the infectivity of the virus to people, and can be affected by the concentration of the virus. They speculated that the infection process follow the nonlinear saturation pattern. They draw the conclusion of epidemic evolution by analyzing the equations in phase space, which provided a reference for public health policies. The validity of the model is verified by comparing it with available data from Bari town in Italy.

In fact, system (1.1) can be modified to simulate other epidemics with oral-faecal transmission ^[6], such as typhoid fever, infectious hepatitis, poliomyelitis, etc. Moreover, system (1.1) can also be used to investigate the spreading of man-environment epidemics ^[6,11,12]. The man-environment epidemics refer that, the infected people may increase the concentration of the virus in the environment, thus people maybe infected by eating the contaminated food ^[11].

Scientific researches show that the virus is highly sensitive to its environment. The concentration of the virus is not evenly distributed in space thus the effect of diffusion cannot be ignored. In 1997, Capasso and Wilson ^[11] studied the following diffusive epidemic model:

where $ u(x, t), \; v(x, t) $ represent the spatial concentration of the bacteria and the infected people, respectively. $ \alpha_1 $, $ \alpha_2 $ are the same as they are in (1.1). For system (1.2), they performed a detailed analysis of the steady-state bifurcation at the endpoints of a one-dimensional interval based on the monotone techniques. They speculated that the system (1.2) has saddle point structure in the natural function space, similar to the ODE case in which the diffusivity of bacteria is also set to zero.

In the transmission of the epidemics, an important factor that can not be ignored is the delay ^[13,14,15]. For example, there is time delay between the release of the virus from the infected person and its successful entry into the susceptible people. For those infected people, it also takes time from being infected to being able to release the virus. All of these can lead to time-delay feedback in the epidemics ^{[16,17,18,19]}. Time delay usually destroys the stability of the solutions, devotes to the complex dynamics of the system. Recently, systems with two delays have been used in characterizing the phenomena in various fields, such as power systems and predator-prey systems ^{[20,21,22,23,24]}, etc.. Among them, Wu and Hsu ^[24] proposed a monostable epidemic model with double delays as follows:

They studied the existence of entire solutions in the system (1.3) with and without the quasi-monotone condition. Some special cases of system (1.3) have been extensively investigated. For example, when $ \tau_1 = d_2 = 0 $, $ \tau_2 $, $ d_1 > 0 $, Zhao and Wang ^[25] studied the properties of the existence and non-existence of the traveling wave fronts; when $ \tau_1 = \tau_2 = d_2 = 0 $, $ d_1 > 0 $, Xu and Zhao ^[26] proved the existence, uniqueness and global exponential stability of the bistable traveling wave fronts; Hsu et al. ^[27] studied the existence and exponential stability of traveling wave solutions in an extended form of system (1.3).

Motivated by the previous works, we study an epidemic model with two delays subjecting to the Neumann boundary conditions in the following form:

where $ u(x, t) $ and $ v(x, t) $ denote the concentration of the virus and the number of the infected people, respectively. $ \frac{1}{\alpha_1} $ denotes the average lifetime of the infectious agent, $ \frac{1}{\alpha_2} $ denotes the average infectious period of infected people, $ S(v) $ represents the contribution of infected people to the growth of the virus, $ T(u) $ represents the infection rate of the virus to people. $ d_1, d_2 > 0 $ are the diffusion coefficients, $ \tau_1 $ represents the average time taken by the infected people from being infected to being able to release the virus, $ \tau_2 $ represents the average time taken by the virus from being released to entrying the susceptible people and $ \overline{\tau} = \max\{\tau_1, \tau_2\} $.

In the spread of an epidemic, the concentration of the virus would increase with the number of the infected people, and the susceptible people would more likely to be infected. Therefore, $ S(v) $ and $ T(u) $ are increasing functions of $ v $ and $ u $ respectively in general. However, the increasing number of the virus poses a threat to public health, people will take effective measures to avoid them from being infected by the virus, including isolating the infected people, regular disinfection, reminding people to pay attention to hygiene through media reports and so on ^[28,29]. Scientific studies have shown that these measures are effective in reducing infection rates ^[30]. Once the virus exceeds a certain concentration $ u_0 $, even if the concentration continues to increase, people begin to take certain measures to reduce the infectivity of the virus. So in this paper, we assume that once if $ u > u_0 $, $ T(u) $ is a decreasing function with respect to $ u $.

In this paper, we study the existence of the positive constant steady state. By analyzing the zeros of characteristic equations of (1.4), we investigate the local stability of constant steady state. We found that when the sum of two delays reaches a critical value, the system undergoes Hopf bifurcation. By taking one of two delays as a bifurcation parameter, we calculate the normal form on the center manifold near the Hopf singularity.

This paper is organized as follows. In section 2, the stability of the positive constant steady state of the system (1.4) and the existence of Hopf bifurcation are established by analyzing the distribution of the eigenvalues. In section 3, the normal form on the center manifold near the Hopf bifurcation point is derived. In section 4, we carry out some numerical simulations to illustrate the theoretic results.

2.   Stability of the positive constant steady state and existence of Hopf bifurcation

In this section, we will study the local stability of the positive constant steady state and the existence of Hopf bifurcation by analyzing the distribution of the eigenvalues.

For $ S(v) $ and $ T(u) $ in (1.4), we make the following assumption:

Normally speaking, the contribution function $ S(v) $ is sufficiently smooth with respect to the infected people $ v $, the infection rate $ T(u) $ is sufficiently smooth with respect to $ u $. For the research purposes, we assume $ S, T\in C^3(\mathbb{R}, \mathbb{R}) $. Particularly, if no one is infected ($ v = 0 $), then no infected people contribute to the growth of the virus, i.e., $ S(0) = 0 $. Similarly, if the concentration of the virus is zero($ u = 0 $), then no one can be infected by the virus, i.e., $ T(0) = 0 $. In the early stage of a transmission, the virus concentration is low but the number of the susceptible people is large, thus the epidemic would spread at a certain speed; as the epidemic continues to spread, although the virus concentration is high, a reduction in the number of the susceptible people will limit the spread of the epidemic. Thus the virus and the population may reach a coexistence state. In $ {\bf (H_1)} $, we denote the number of infected people in this coexistence state as $ v^* $.

Define $ u^* = S(v^*)/\alpha_1 $, under the assumption $ (\bf H_1) $ system (4) has two constant steady states $ E_0(0, 0) $ and $ E_*(u^*, v^*) $. Obviously, $ E_* $ is a positive constant steady state. To investigate the dynamics of the epidemics, we focus on the properties of $ E_* $.

the linearized equation of system (1.4) at $ E_* $ is

where

and $ u(x, t) $, $ v(x, t) $ still satisfy the homogeneous Neumann boundary conditions.

The characteristic equations of the system (2.1) are given by

where $ I_2 $ is a 2 $ \times $ 2 identity matrix, $ M_n = -\frac{n^2}{l^2}{D} $. Since

denote $ m = S_{v^*}\times T_{u^*} $, Eq (2.2) becomes

Denote $ \tau = \tau_1+\tau_2 $, then Eq (2.3) becomes

For convenience, denote $ k_1 = d_1\frac{n^2}{l^2}+\alpha_1, \; k_2 = d_2\frac{n^2}{l^2}+\alpha_2 $, then we have $ k_1, k_2 > 0 $. Now the characteristic Eq (2.4) with $ \tau = 0 $ is in the following form:

To ensure all the roots of Eq (2.5) have negative real parts, we further make the following assumption

Then we can conclude the following lemma.

Lemma 2.1. When $ \tau = 0 $, all the roots of Eq (2.4) have negative real parts under the assumption $ {\bf(H_2)} $.

Since the proof of Lemma 2.1 is intuitive, here we omit it.

Now we will study the stability of $ E_* $ and investigate the existence of Hopf bifurcation of system (1.4). Assume that $ \lambda = {\rm i}\omega\; (\omega > 0) $ is a root of (2.4). Substituting $ {\rm i}\omega $ into (2.4) and separating the real and imaginary parts leads to

This yields

Denote $ z = \omega^2 $, then the above equation becomes

the roots of (2.7) can be given by

Clearly, $ z_2 < 0 $. To maintain the existence of a positive $ \omega $ that satisfies (2.6), $ z_1 > 0 $ should hold. Obviously, $ z_1 > 0 $ provided that $ m^2 > k_1^2k_2^2 $.

If $ m^2\leq \alpha_1^2 \alpha_2^2 $, then for all $ n\in\mathbb N_0 $, we have $ m^2\leq (d_1\frac{n^2}{l^2}+\alpha_1)^2(d_2\frac{n^2}{l^2}+\alpha_2)^2 $, and thus $ z_1\leq 0 $, all roots of (2.4) have negative real parts. Moreover, if $ m < \alpha_1\alpha_2 $ and $ m^2\leq \alpha_1^2 \alpha_2^2 $ are satisfied, i.e., $ -\alpha_1\alpha_2 \leq m < \alpha_1\alpha_2 $, then $ E_* $ is locally asymptotically stable. We can draw the following proposition.

Proposition 1. If $ -\alpha_1\alpha_2 \leq m < \alpha_1\alpha_2 $, all roots of (2.4) have negative real parts, then $ E_* $ is locally asymptotically stable.

Now we further make the following assumption:

Then we study the existence of purely imaginary roots of (2.4) in the following lemma.

Lemma 2.2. If $ {\bf(H_3)} $ is satisfied, then Eq (2.4) with $ n\in\{0, 1, \cdots, N-1\} $ has a pair of purely imaginary roots $ \pm {\rm i}\omega_n $ when $ \tau = \tau_{j, n}, $ where

Proof. From $ {\bf(H_3)} $ we know that Eq (2.7) with $ n\in\{0, 1, \cdots, N-1\} $ has a positive root denoted by $ z_{n} $. Hence, $ \omega_n = \sqrt{z_{n}} > 0 $ makes sense. Define

From the second equation of (2.6) one obtains that $ \tau_{0, n}\omega_n\in(0, \pi] $. Then $ (\omega_n, \tau_{j.n}) $ solves (2.6), which implies that $ \pm{\rm i}\omega_n $ is a pair of purely imaginary roots of Eq (2.4) with $ \tau = \tau_{j, n} $.

Lemma 2.3. If assumptions $ {\bf(H_2)} $ and $ {\bf (H_3)} $ are satisfied, then $ \omega_n $ is monotone decreasing with respect to $ n $.

Proof. By Lemma 2.2, we know that $ \omega_n $ exists for $ n = 0, \; 1, \cdot\cdot\cdot\; N-1 $. Assumptions $ {\bf(H_2)} $ and $ {\bf(H_3)} $ imply that $ m < 0 $. We know that

Hence, $ z_n $ is monotone decreasing with respect to $ n $. For $ i = 0, 1, \cdots, N-1 $ we have $ z_{i} > z_{i+1} $. Thus $ \omega_n $ is monotone decreasing with respect to $ n $.

Lemma 2.4. If assumption $ {\bf(H_2)} $ and $ {\bf(H_3)} $ are satisfied, then $ \tau_{0, n}\omega_n $ is monotone increasing with respect to $ n $. Moreover, $ \tau_{0, 0} = \min\{\tau_{0, n}\}_{0\leq n\leq N-1} $.

Proof. From Lemma 2.3, for $ n = 0, 1, \cdots, N-2 $, we have $ \omega_n > \omega_{n+1} $ and

By $ -{\omega_n}^2 < -{\omega^2_{n+1}} $ and $ m < 0 $, we have $ \dfrac{-{\omega_n}^2+(d_1\frac{n^2}{l^2}+\alpha_1)(d_2\frac{n^2}{l^2}+\alpha_2)}{m} $ is monotone decreasing with respect to $ n $. Since the anti-trigonometric function $ y = \arccos x $ is monotone decreasing with respect to $ x $, from the properties of composite functions, we can obtain that $ \tau_{0, n}\omega_n $ is monotone increasing with respect to $ n $. The last conclusion follows from $ \omega_n > \omega_{n+1} $.

Under the assumption $ {\bf (H_2)} $ and $ {\bf (H_3)} $, let $ \lambda(\tau) = \alpha(\tau)+{\rm i}\beta(\tau) $ be the root of Eq (2.4) with $ n\in\{0, 1, \cdots, N-1\} $ satisfying $ \alpha(\tau_{j, n}) = 0 $ and $ \beta(\tau_{j, n}) = \omega_n $.

Lemma 2.5. Suppose that $ {\bf (H_2)} $ and $ {\bf (H_3)} $ are satisfied. Then

Proof. Since Eq (2.4) can be rewritten in the following form:

Substituting $ \lambda(\tau) $ into (2.9) and taking the derivative of both sides of Eq (2.9) with respect to $ \tau $ gives

Then we have

Hence,

This completes the proof. By Lemmas 2.1, 2.4, 2.5 and applying Corollary 2.4 in ^[31], we have the following conclusions on the distribution of the roots of the characteristic Eq (2.4).

Lemma 2.6. Suppose that $ {\bf(H_2)} $ and $ {\bf(H_3)} $ are satisfied. Then there exists a sequence values of $ \left\{\tau_{j, n}\right\}(j = 0, 1, 2, \cdots; n = 0, 1, \cdots, N-1) $, such that all the roots of Eq (2.4) have negative real parts when $ \tau_1+\tau_2\in[0, \tau_{0, 0}) $; When $ \tau_1+\tau_2 = \tau_{0, 0} $, all the other roots of Eq (2.4) have negative real parts, except for a pair of pure imaginary roots $ \pm{\rm i}\omega_0 $; When $ \tau_1+\tau_2 > \tau_{0, 0} $, Eq (2.4) has at least a couple of roots with positive real parts.

From Lemma 2.6, we have the following the conclusions on the stability of the positive equilibrium of the system (1.4) and the existence of Hopf bifurcations directly.

Theorem 2.7. Suppose that $ {\bf(H_1), (H_2)} $ and $ {\bf(H_3)} $ are satisfied. Then the positive constant steady state $ E^* $ of system (1.4) is asymptotically stable when $ \tau_1+\tau_2\in[0, \tau_{0, 0}) $, unstable when $ \tau_1+\tau_2 > \tau_{0, 0} $, and system (1.4) undergoes a Hopf bifurcation at $ E_* $ when $ \tau_1+\tau_2 = \tau_{j, n} $, $ j = 0, 1, 2, \cdots; \; \; n = 0, 1, \cdots, N-1. $

3.   Normal form of the Hopf bifurcation on the center manifold

In this section, we will calculate the normal form of Hopf bifurcation on the center manifold. To investigate the dynamical behaviours of system (1.4) near the Hopf bifurcation point, we will study the normal form of Hopf bifurcation by making use of the normal form method of partial functional differential equations ^[32,33].

Without loss of generality, we assume that $ \tau_1 > \tau_2 $, let $ \overline u(x, t) = u(x, \tau_1t)-u^* $, $ \overline v(x, t) = v(x, \tau_1t)-v^* $ and take out of the bars, the system (1.4) becomes

with the boundary and initial conditions:

where

We define the real-valued Hilbert space

and the complexification of $ Y $

with the general complex-value $ L^2 $ inner product

for $ U = (u_1, u_2)^T, V = (v_1, v_2)^T \in Y_\mathbb{C} $.

Let $ \mathcal{C}: = C([-1, 0], Y_\mathbb{C}) $ represents the phase space with the sup norm. We write $ u_t\in\mathcal{C} $ for $ u_t(\theta) = u(t+\theta) $, $ -1\leq\theta\leq0 $.

In the previous section, we have obtained that when $ \tau = \tau_{0, 0} $, all the roots of Eq (2.4) have negative real parts except $ \pm{\rm i}\omega_0 $. We denote the Hopf bifurcation point by $ (\tau_1^*, \tau_2^*) $ with $ \tau_1^*+\tau_2^* = \tau_{0, 0} $, and introduce parameter $ \mu $ by setting $ \tau_1 = \tau_1^*+\mu $ with fixing $ \tau_2 = \tau_2^* $ and $ |\mu| $ sufficiently small. Then we have $ \tau = \tau_1+\tau_2 = \tau_1^*+\tau_2^*+\mu = \tau_{0, 0}+\mu $, and $ \mu = 0 $ is the Hopf bifurcation value of system (3.1).

Denote $ U(t) = (u(x, t), v(x, t))^T $ and substitute $ \tau_1 = \tau_1^*+\mu, \tau_2 = \tau_2^* $ into (3.1). We have

where

and

Consider the linearized system of (3.3)

In $ Y $, the eigenvalues of $ D\Delta $ are $ -d_1\frac{n^2}{l^2} $ and $ -d_2\frac{n^2}{l^2} $, $ n\in\mathbb N_0 $, with the corresponding normalized eigenvectors: $ \beta_n^{(1)}(x) = \gamma_n(x)e_1 $, $ \beta_n^{(2)}(x) = \gamma_n(x)e_2 $, where

and $ e_1 = (1, 0)^T, e_2 = (0, 1)^T $ are the unit coordinate vectors of $ \mathbb R^2 $.

Define the subspace of $ \mathcal{C} $ as $ \mathcal{B}_n $, by

we write $ \left(v(\cdot), \beta_n\right) = \left(\left(v(\cdot), \beta_n^{(1)}\right), \left(v(\cdot), \beta_n^{(2)}\right)\right)^T $ for simplification. Then on $ \mathcal{B}_n $, the linear Eq (3.4) is equivalent to retarded functional differential equation on $ \mathbb{R}^2 $:

By Riese representation theorem, there exists a $ 2\times 2 $ matrix function $ \eta_n(\mu, \theta) $, $ \theta\in [-1, 0] $, whose entries are of bounded variation such that

In fact, we can choose

For simplicity, we write $ \eta_n(\theta) $ for $ \eta_n(0, \theta) $, define $ \mathcal{A} $ as the infinitesimal generator of the semigroup generated by Eq (3.5) with $ \mu = 0 $ and $ n = 0 $. Define the adjoint operator $ \mathcal{A}^* $ of $ \mathcal{A} $ on $ \mathcal{C}^*: = C([0, 1], (\mathbb{R}^2)^*) $ as

with the bilinear form

Let

be the eigenvectors of $ \mathcal{A} $ and $ \mathcal{A}^* $ corresponding to $ {\rm i}\omega_0\tau_1^* $ and $ -{\rm i}\omega_0\tau_1^* $ respectively. Thus

where

We choose $ M $ as

which assures that $ \langle q^*(s), q(\theta)\rangle = 1 $.

Let

and define $ \mathcal{A}_\mu $ as the infinitesimal generator of the $ C_0 $ semi-group of solution maps of the linear equation of (3.3):

with

where

On $ \mathcal{BC} $, Eq (3.3) can be rewritten as an abstract ordinary differential equation

Using the same notations in Wu ^[34], we compute the coordinates to describe the center manifold $ \mathcal C_0 $ at $ \mu = 0 $. Let $ U_t $ be the solutions of Eq (3.3) with $ \mu = 0 $. Define

On the center manifold $ \mathcal C_0 $, we have

where

$ z $ and $ \bar z $ are local coordinates for center manifold $ \mathcal C_0 $ in the direction of $ q^* $ and $ \overline q^* $ respectively. Notice that $ W $ is real if $ U_t $ is real, we consider real solutions.

For the solution $ U_t\in\mathcal C_0 $ of the system (3.3) with $ \mu = 0 $,

It can be denoted as

where

Base on (3.7) and (3.8), similar to the method used in Zhao and Wei ^[35], we obtain

where

By expanding the above series and comparing the coefficients, we have

Notice that

where

and

We have

Comparing the coefficients with (3.9), we obtain that

We need to compute $ W_{20}(\theta) $ and $ W_{11}(\theta) $ for $ \theta\in[-1, 0) $, since

Comparing the coefficients with (3.9), we obtain

From (3.10), it can be given that

Solving for $ W_{20}(\theta) $, we have

and similarly

where $ E_1 $ and $ E_2 $ can be determined by setting $ \theta = 0 $ in $ H $. Since

we have

From (3.10) and the definition of $ \mathcal A_0 $, we have

Noting that

substituting (3.12) into (3.16), we obtain

Substituting (3.14) into this, we have

Solving Eq (3.17) for $ (E_1^{(1)}, E_2^{(2)})^T = E_1 $, we get

Similarly, we can get

and hence,

Then $ g_{21} $ can be confirmed. Therefore, we can calculate the following quantities

By the general theory in Hassard ^[36], we obtain that $ \mu_2 $ determines the direction of the Hopf bifurcation: When $ \mu_2 > 0 (<0) $, the direction of the Hopf bifurcation is forward (backward), in other words, for $ \tau_1 > \tau_1^* (\tau_1 < \tau_1^*) $, there exist the bifurcating periodic solutions; $ \beta_2 $ determines that whether the bifurcating periodic solutions are stable: When $ \beta_2 < 0 (>0) $, the bifurcating periodic solutions are orbitally asymptotically stable (unstable); $ T_2 $ determines the period of the bifurcating periodic solutions: When $ T_2 > 0 (<0) $, the period of the bifurcating periodic solutions increases (decreases).

4.   Numerical simulations

In this section, we choose

and take the following data:

Then system (1.4) with (4.1) and (A) has a unique positive constant steady state $ E_* = (16.8, 0.336) $. By direct calculation, we have $ S_{v^*} = 40 $, $ T_{u^*} = -0.068 $ and $ m = -2.72 $. Thus the assumptions $ (\bf{H_1}) $ and $ (\bf{H_2}) $ are satisfied, and $ E_* $ is asymptotically stable when $ \tau_1 = \tau_2 = 0 $.

By calculation, we get that when $ N = 3, \; m^2 > (d_1\frac{n^2}{l^2}+\alpha_1)^2(d_2\frac{n^2}{l^2}+\alpha_2)^2 $ for $ 0\leq n < N\; and\; m^2\leq (d_1\frac{N^2}{l^2}+\alpha_1)^2(d_2\frac{N^2}{l^2}+\alpha_2)^2, $ thus assumption $ (\bf{H_3}) $ is satisfied. From Eq (2.8), we get that

And Eq (2.7) has no positive roots for $ n\geq3 $. Clearly, $ \tau_{0, 0}\approx 0.7023. $

From Theorem 2.7, we know that the positive constant steady state $ E_* $ is asymptotically stable when $ \tau_1+\tau_2\in [0, \tau_{0, 0}) $, and unstable when $ \tau_1+\tau_2 > \tau_{0, 0} $. Meanwhile, (1.4) undergoes a Hopf bifurcation at $ E_* $ when $ \tau_1+\tau_2 = \tau_{j, n} $, $ n = 0, 1, 2; j\in\mathbb N_0 $. The bifurcation set on $ (\tau_1, \tau_2) $ plane is shown by Figure 1.

When $ (\tau_1, \tau_2) $ are chosen as $ P_1(0.4, 0.2) $, we have $ \tau_1+\tau_2 < 0.7023 $, thus $ E_* $ is asymptotically stable, this is illustrated in Figure 2.

We choose the point $ P_3(0.55115, 0.15115) $ on the curve $ \tau = \tau_{0, 0} $, that is $ \tau_1^* = 0.55115 $ and $ \tau_2^* = 0.15115 $. According to the algorithm given in the previous section, we can obtain that $ c_1(0) = -0.00155-0.00290{\rm i} $, $ \lambda'(\tau_1^*) = 0.92063-0.94716{\rm i} $, $ \mu_2 = 0.00168 > 0 $, $ \beta_2 = -0.0031 < 0 $ and $ T_2 = 0.00565 > 0 $. This implies that the direction of the Hopf bifurcation is forward, the bifurcating periodic solutions are stable and the period of the bifurcating periodic solutions increases, which are illustrated by Figure 3.

In fact, with the help of the norm form calculated in section 2, we get that when $ \tau_1^* $ and $ \tau_2^* $ satisfying $ \tau_1^*+\tau_2^* = \tau_{0, 0} $ and $ \tau_1^* > \tau_2^* $, the corresponding $ \text{Re}c_1(0) $ would be negative, we list some values of $ c_1(0) $ in Table 1.

5.   Conclusions

In this paper, we study the dynamics of epidemics by using a diffusive model with two delays. First we perform the stability analysis at the constant steady state in the system. Then we present the conditions for the existence of the purely imaginary roots of the characteristic equations. Moreover, we investigate the existence of Hopf bifurcation by taking two delays as parameters. Theoretical result shows that the sum of two delays affect the stability of the steady state and the existence of Hopf bifurcation. Considering the sum of two delays as one parameter, we verify that $ \tau_1^*+\tau_2^* = \tau_{0, 0} $ is the critical Hopf bifurcation curve and the stability switching curve of this system. To calculate the normal form on the center manifold near the Hopf singularity, we fix $ \tau_2 $ as a constant and take $ \tau_1 $ as a parameter. Finally, numerical simulations are carried out to verify the theoretical results. The simulations show that when $ \tau_1+\tau_2 < \tau_{0, 0} $, $ E_* $ is locally asymptotically stable; when $ \tau_1+\tau_2 > \tau_{0, 0} $, the bifurcating periodic solutions are stable. Both theoretical verifications and numerical simulations reveal that when $ \tau_1+\tau_2 < \tau_{0, 0} $, the system is locally stable at $ E_* $; when $ \tau_1+\tau_2 > \tau_{0, 0} $ and close enough to $ \tau_{0, 0} $, $ E_* $ is unstable, the periodic solutions generated by the Hopf bifurcation is stable. All the analysis shows that the sum of two delays plays an important role in the dynamics of the system.

Acknowledgments

The authors are grateful to the anonymous referees for their helpful comments and valuable suggestions which have improved the presentation of the paper. This research is supported by National Natural Science Foundation of China (No.11771109).

Conflict of interest

The authors declared that they have no conflicts of interest to this work. We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.