Bifurcation analysis and optimal control of SEIR epidemic model with saturated treatment function on the network

1.   Introduction

Many practical problems in the real world can be abstracted into complex network models for research. For example, the spread and control of epidemics, the spread of computer viruses in the Internet, and the spread of rumors can all be regarded as spreading behaviors that obey a certain law on a complex network. At present, there have been many research results on the spread and control of epidemics on complex networks. For related research reports, see references ^{[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]}.

In classical epidemic-disease dynamics models, it is generally assumed that the rate of treatment for a disease is proportional to the number of infected persons. This means that medical resources such as drugs, vaccines, hospital beds and isolation facilities are sufficient for the epidemic. In reality, however, every community or country has adequate or limited capacity for treatment and vaccination. If too much medical resources are invested, social resources will be wasted. If fewer resources are put into care, the risk of disease outbreaks increases. Therefore, it is important to identify different capacity for treatment depending on the community or country.

In order to investigate the effect of the limited capacity for treatment on the spread of infectious disease, references ^[24,25] established an infectious disease model under a discontinuous treatment strategy. Wang et al. considered the following segments of treatment functions in the literature ^[25]

where $k$ is the cure rate and $I_{0}$ is the maximum capacity of the medical system. In other words, when the number of patients is small and does not exceed the maximum capacity of the medical system, the treatment rate is proportional to the number of patients.When the number of patients exceeds the maximum capacity of the medical system, the treatment rate of the disease is a constant constant. reference ^[26], Zhang and Liu introduced the following saturation processing functions:

where $r$ is the cure rate, $\alpha$ is used to describe the impact of delayed treatment of a person with an illness in a situation where medical care is limited. Related studies on the function of saturation therapy are reported in references ^{[27,28,29,30]}.

This paper proposes a SEIR model with a saturated processing function based on a complex network. We will analyze the influence of the heterogeneity of the population on the spread of the virus, and use the advancement of optimal control theory to study the time-varying control of infectious diseases.

The organization structure of this article is as follows. In Section 2, we propose a network-based SEIR model with saturation processing function. In Section 3, we analyzed the dynamics of the model. In Section 4, we study the bifurcation behavior of the model. In Section 5, we study the optimal control problem. Finally, in Section 6 and Section 7, we give the results of the numerical simulation and summarize the paper.

2.   The model

In this section, we propose a network-based SEIR epidemic model with saturation treatment function. Then we use the SEIR epidemic model to prove the positivity and boundedness of understanding.

In the classic compartmental model, each individual has the same probability of contact with an infected individual. However, when the population is large, it is generally believed that factors such as exposure heterogeneity should be considered. Therefore, a complicated network is added to the infectious disease model to describe contact and so on. In a complex network work, each individual corresponds to a node of the network, and the interaction between two individuals corresponds to a link between two nodes.

On the basis of reference ^[31], we consider the following network-based SEIR epidemic model.

The flow diagram of the state transition is depicted in Figure 1.

The dynamic equations can be written as

The explanation of the parameters is the same as in reference ^[31]. $\theta(t)$ describes a link pointing to an infected node, which satisfies the relation $\theta(t) = \frac{1}{\langle k \rangle}\sum\limits_{k}kp(k)I_{k}(t)$, $p(k)$ is a degree distribution, $\langle k\rangle = \sum\limits_{k = 1}^nkp(k)$ describes the average degree. Here, we assume that the connectivities of nodes in network at each time are uncorrelated.

$\frac{rI_{k}(t)}{1+\alpha\theta(t)}$ represents the recovery of the $k$-th infection group after treatment.

$g(\theta) = \frac{1}{\langle k \rangle}\sum\limits_{h = 1}^nhP(h)\frac{rI_{h}(t)}{1+\alpha\theta(t)} = \frac{r\theta(t)}{1+\alpha\theta(t)}$ represents the probability that a given edge is connected to an infected node that is recovering through treatment. Similar to the processing function given in reference ^[26], g is a saturation processing function.

In the following cases, we make the following assumptions: (ⅰ) all new births are susceptible; (ⅱ) the total number of nodes is constant, and the number of deaths is equal to the number of births, so $\pi = \mu$; (ⅲ) assume that the degree of each node is time-invariant.

$N_{k}$ is a constant, stands for the number of nodes with degree k. Then $N_{k} = S_{k}+E_{k}+I_{k}+R_{k}, k = 1, 2, · · ·, n$ and $\sum\limits_{k}N_{k} = N$.

For the practice, the initial condition for model (2.1) satisfy:

The probability $0 \leq\theta(t)\leq1$ describes a link pointing to an infected host, which satisfies the relation

At t time, the global ensity is

Because

For the sake of convenience. it is assumed that the system has been normalized, so there are $S_{k}(t)+E_{k}(t)+I_{k}(t)+R_{k}(t) = 1$, So the positive invariant set of system (3.1) is

3.   Global behavior of the model

3.1.   Equilibria and basic reproduction number

Lemma 1. Suppose that $S_k(t), E_k(t), I_k(t), R_k(t)$ is a solution of model $(2.1)$ satisfying initial conditions of Eq.$(2.2)$. Then $\Omega = \{(S_{1}, E_{1}, I_{1}, R_{1}\cdots, S_{n}, E_{n},$ $I_{n}, R_{n}):S_{k}\geq0, E_{k}\geq0,$ $I_{k}\geq0, R_{k}\geq0, S_{k}+E_{k}+I_{k}+R_{k} = 1, k = 1, 2, \cdots, n\}$ is a positively invariant for model $(2.1)$.

Proof. First, we proved that $E_{k}(t)\geq 0$ for any $t > 0$. If not, as $E_{k}(0) > 0$, there exist $k_{0}\in\{1, 2, · · ·, n\}$ and $t^{\ast}$ such that $t^{\ast} = \inf\{t\arrowvert E_{k_{0}}(t) = 0, \dot{E}_{k_{0}}(t) < 0\}$. On the hand, according to the definition of $t^{\ast}$, we have $E_{k_{0}}(t^{\ast}) = 0, \dot{E}_{k_{0}}(t^{\ast}) < 0$ and $E_{k_{0}}(t) > 0$ for $0\leq t < t^{\ast}$. From the second equation of model (2.1), we obtain $\dot{E}_{k_{0}}(t^{\ast}) = \beta k_{0}\lambda S_{k_{0}}(t^{\ast})\theta(t^{\ast}) < 0$, which leads to $S_{k_{0}}(t^{\ast}) < 0$. On the other hand, in the time interval $[0, t^{\ast}]$, from the first equation of model (3.1), we have $\dot{S}_{k_{0}}(t)\geq -\beta kS_{k_{0}}(t)\theta(t)-\mu S_{k_{0}}(t)$, which implies $S_{k_{0}}(t)\geq S_{k_{0}}(0) e^{-\beta k\theta(t)-\mu}\geq0$ for $0\leq t\leq t^{\ast}$. In particular, it follows $S_{k_{0}}(t^{\ast})\geq0$. It is a contradiction, so $E_{k}(t)\geq 0$ for any $t > 0$. Using the same method, we can verify $S_{k}(t)\geq0$, $I_{k}(t) \geq0$ and $R_{k}(t) \geq0$ for any $t > 0$.

The proof is complete.

Theorem 1. Basic reproduction number of system $(2.1)$ is $R_{0} = \frac{\beta\delta}{(\mu+\delta)(r+\mu+\gamma)}\frac{\langle k^2\rangle}{\langle k\rangle}$. If $R_{0} < 1$, the system $(2.1)$ has only virus-free quilibrium; if $R_{0} > 1$, then system $(2.1)$ has endemic equilibrium.

Proof. Assuming that $E = (S_{1}, E_{1}, I_{1}, R_{1}, \cdots, S_{n}, E_{n}, I_{n}, R_{n})$ is the unique balance of system (2.1), then E should satisfy

The equilibria of system $(2.1)$ are determined by setting

After calculation, there is

This implies that

Substituting $I_{k}$ into system $(2.3)$, we have

$\theta[1-\frac{1}{\langle k \rangle}\sum\limits_{h = 1}^n\frac{\beta h^{2}P(h)}{\frac{p\beta h\theta+\mu+\delta}{p\beta h\theta+\delta}(\frac{r}{1+\alpha\theta}+\mu+\gamma)+\frac{\beta h\theta}{p\beta h\Theta+\delta}(\frac{r}{1+\alpha\theta}+\mu+\gamma)+\beta h\theta+\frac{\beta h\theta}{\mu}(\frac{r}{1+\alpha\theta}+\gamma)}] = 0$

Define

$F(\theta) = 1-\frac{1}{\langle k \rangle}\sum\limits_{h = 1}^n\frac{\beta h^{2}P(h)}{\frac{p\beta h\Theta+\mu+\delta}{p\beta h\theta+\delta}(\frac{r}{1+\alpha\theta}+\mu+\gamma)+\frac{\beta h\theta}{p\beta h\theta+\delta}(\frac{r}{1+\alpha\theta}+\mu+\gamma)+\beta h\theta+\frac{\beta h\theta}{\mu}(\frac{r}{1+\alpha\theta}+\gamma)}$

Because

and $F$ is continuous on $[0, 1]$. According to the intermediate value theorem, if $F(0) < 0$, then $F(\theta) = 0$ has a positive solution. While

namely

Then let $R_{0} = \frac{\beta\delta}{(\mu+\delta)(r+\mu+\gamma)}\frac{\langle k^2\rangle}{\langle k\rangle}$. If $R_{0} < 1$, the system $(2.1)$ has only disease-free quilibrium $E_{0}$; if $R_{0} > 1$, then system (1) has endemic equilibrium.

3.2.   Stability of disease-free equilibrium

Remark 1. One can also obtain the same basic reproduction number by using the method of second generation matrix ^[32], where indicates that if $R_{0} < 1$, then the disease-free equilibrium $E_{0}$ of model $(2.1)$ is locally asymptotically stable; if $R_{0} > 1$, then $E_{0}$ is unstable.

Theorem 2. For system $(2.1)$, if $R_{0} < 1$, then the free quilibrium $E_{0}$ is locally asymptotically stable.

Proof. $1.$ Since $S_{k}+E_{k}+I_{k}+R_{k} = 1$ system $(2.1)$ is converted into an equivalent system:

The Jacobian matrix of model $(3.5)$ at disease-free equilibrium $J_{2n\times2n}$ is given by:

and

The characteristic equation of the disease-free equilibrium is

The eigenvalues of $J_{E_{0}}$ are all negative when $R_{0} < 1$, the disease-free equilibrium $E_{0}$ of model $(2.1)$ is locally asymptotically stable when $R_{0} < 1$.

Theorem 3. For system $(2.1)$, Denote $\hat{R}_{0} = \frac{\beta\langle k^{2}\rangle} {(\frac{r}{1+\alpha}+\mu+\gamma)\langle k\rangle}$. If $\hat{R}_{0} < 1$, then the free quilibrium $E_{0}$ is globally asymptotically stable.

Proof. Constructing a lyapunov function:

Deriving the V function along system $(2.1)$,

We now claim that if $0 < S_{k}(t) < 1$, the first term is negative. In fact, because of $\pi = \mu$, we have

This means

Because $0 < I_{k}(t) < 1$ for all $k$, we know that $0 < \theta(t) < 1$. Therefore

If $\hat{R}_{0} < 1$, $\left.\frac{dV}{dt}\right|_{(2.1)}\leq 0$. The equation holds if and only if $I_{k}(t) = 0$. For the limit system $E^\prime_{k}(t) = -(\mu+\delta)E_{k}(t)$, it is easy to see that $\lim\limits_{t\rightarrow +\infty}E_{k}(t) = 0$, For the limit system $R^\prime_{k}(t) = -\mu R_{k}(t)$, it is easy to see that $\lim\limits_{t\rightarrow +\infty}R_{k}(t) = 0$.Finally, from $S_{k}(t)+E_{k}(t)+I_{k}(t)+R_{k}(t) = 1$, we get $\lim\limits_{t\rightarrow +\infty}S_{k}(t) = 1$. According to the LaSalle invariant principle, $E_{0}$ is globally attractive. Combining the locally asymptotically stability, we conclued that the free quilibrium $E_{0}$ is globally asymptotically stable.

We noticed that in Theorem 3, the disease-free equilibrium $E_{0}$ is globally asymptotically stable only when $\hat{R}_{0} < 1$. In addition, you can see that $R_{0} < \hat{R}_{0}$. This prompts us to think about whether the disease can persist even if $R_{0} < 1$.

4.   Bifurcation analysis

In this section, we will determine whether there is a backward bifurcation at $R_{0} = 1$. More accurately, we derive the condition for determining the bifurcation direction at $R_{0} = 1$. let us find endemic equilibria. At an endemic equilibrium, we have $I_{k} > 0$ for all $k = 1, 2, ..., n$, which means $\theta > 0$. From $F(\theta) = 0$, we can see that the unique equilibrium should satisfy the following equation :

The denominator and numerator are multiplied by the quantity $\frac{\delta\langle k^{2}\rangle}{(\mu+\delta)(r+\mu+\gamma+\mu_{d})\langle k\rangle}$, which can be expressed by $R_{0}$ and $\theta$ as follows:

$\frac{1}{\langle k \rangle}\sum\limits_{h = 1}^n\frac{R_{0}h^{2}P(h)}{\frac{\delta\langle k^{2}\rangle}{(\mu+\delta)(r+\mu+\gamma)\langle k\rangle}\frac{p\beta h\theta+\mu+\delta}{p\beta h\theta+\delta}(\frac{r}{1+\alpha\theta}+\mu+\gamma)+\frac{R_{0} h\theta}{p\beta h\theta+\delta}(\frac{r}{1+\alpha\theta}+\mu+\gamma+\mu_{d})+R_{0} h\theta+\frac{R_{0} h\theta}{\mu}(\frac{r}{1+\alpha\theta}+\gamma)} = 1$

If the local $\theta$ is a function of $R_{0}$, the sign of the bifurcation direction is the slope at $(R_{0}, \theta) = (1, 0)$ (cf. Figure 2). More specifically, if the derivative is positive at the critical value $(R_{0}, \theta) = (1, 0)$, that is

the endemic equilibrium curve diverges forward. Conversely, if the derivative is negative at the critical value $(R_{0}, \theta) = (1, 0)$, that is

the local equilibrium curve diverges backward.

$\frac{1}{\langle k \rangle}\sum\limits_{h = 1}^n\frac{I_{1}-I_{2}}{\{\frac{\delta\langle k^{2}\rangle}{(\mu+\delta)(r+\mu+\gamma)\langle k\rangle}\frac{p\beta h\theta+\mu+\delta}{p\beta h\theta+\delta}(\frac{r}{1+\alpha\theta}+\mu+\gamma)+\frac{R_{0} h\theta}{p\beta h\theta+\delta}(\frac{r}{1+\alpha\theta}+\mu+\gamma)+R_{0} h\theta+\frac{R_{0} h\theta}{\mu}(\frac{r}{1+\alpha\theta}+\gamma)\}^{2}} = 0$

where

where

Substituting $R_{0} = 1$ amd $\theta = 0$ into Eq (4.1), we have

This yields

From this we have

We can from Eq (4.2) that

To conclude, we have the following theorem:

Theorem 4. System (2.1) has backward bifurcation at $R_{0} = 1$ if and only if

where $\langle k^{2}\rangle = \sum\limits_{h = 1}^nh^{2}P(h)$, $\langle k^{3}\rangle = \sum\limits_{h = 1}^nh^{3}P(h)$.

It can be seen from Theorem 4 that the nonlinear processing function does play a key role in causing backward bifurcation. More precisely, if $\alpha$ is large enough to satisfy condition (4.3), backward bifurcation will occur. Otherwise, when this effect is weak, there is no backward bifurcation.

5.   Optimal quarantine control

In order to achieve the control goal and reduce the control cost, the optimal control theory is a feasible method.

Due to practical needs, we define a limited terminal time. Then, model $(2.1)$ is rewritten as

The objective function is set as

with Lagrangian

Our objective is to find a optimal control $r_{k}^\ast(t)$ such that

where $U = \{r_{k}, 0\leq r_{k}(t)\leq 1, t\in [0, T]\}$ is the control set. Next we will analyze the optimal control problem.

5.1.   The existence of optimal control solution

In order to obtain the analytic solution of the optimal control, we need to prove the existence of the optimal control first. The following lemma comes from reference ^[33].

Lemma 2. If the following five conditions can be met simultaneously:

(C1) $U$ is closed and convex.

(C2) There is $r\in U$such that the constraint dynamical system is solvable.

(C3) $F(X(t), r_{k}(t))$ is bounded by a linear function in $X$.

(C4) $L(X(t), r_{k}(t))$ is concave on $U$.

(C5) $L(X(t), r_{k}(t))\geq c_{1}\left\|r\right\|_2^\varphi+C_2$ for some $\varphi > 1$, $c_{1} > 0$ and $c_{2}$.

Then the optimal control problem has an optimal solution.

Theorem 5. There exists an optimal solution $r^* (t)$ stastify the control system.

Proof. Let us show that the five conditions in Lemma 3 hold true.

ⅰ) It is easy to prove the control set $U$ is closed and convex. Suppose $r$ is a limit point of $U$, there exists a sequence of points $\{r_{n}\}_{n = 1}^\infty$, we have $r\in (L^2[0, T])^n$. The closedness of $U$ follows from $0\leq r = \lim\limits_{n\rightarrow \infty}r_{n}\leq 1$.

Let $r_1, r_2\in U$, $\eta\in (0, 1)$, we have $0\leq (1-\eta)r_1+\eta r_2\leq 1\in (L^2[0, T])^n$ as $(L^2[0, T])^n$ is a real vector space.

ⅱ) For any control variable $r\in U$, the solution of system $(12)$ obviously exists following from the Continuation Theorm for Differential Systems ^[34].

ⅲ) Let the function $F(X(t), r_{k}(t))$ resprensts the right side of model $(5.1)$, $F$ is continuous, bounded and can be written as a linear function of $X$ in three state.

ⅳ) $L(X(t), r_{k}(t))$ is concave on $U$. We calculate the function $\frac{\partial^2 L}{\partial r^2} = r_{k}(t)\geq 0$, so the function L is concave on $U$.

ⅴ) Since $I_{k}(t)\geq 0$, $L\geq\sum\limits_{k = 1}^n\frac{1}{2}A_{k}r^2_{k}(t)\geq\frac{1}{2}min\{A_{k}\}\left\|r\right\|_2^2$. There exsits $\varphi = 2$, $C_1 = \frac{1}{2}min\{A_{k}\}$ and $C_2 = 0$ such that $L(X(t), r_{k}(t))\geq c_{1}\left\|r\right\|_2^\varphi+C_2$.

5.2.   Solution to the optimal control problem

In this section, we will deal with the OCP based on the Pontryagin Maximum Principle ^[35]. Define the Hamiltonian $H$ as

where $\lambda_{1k}(t), \lambda_{2k}(t), \lambda_{3k}(t), \lambda_{4k}(t)$ are the adjoint variables to be determined later.

Theorem 6. Let $S^*_{k}(t)$, $E^*_{k}(t)$, $I^*_{k}(t)$, $R^*_{k}(t), k = 1, 2, \cdots, n$ be the optimal state solutions of dynamic model $(5.1)$ related to the optimal control $r^* (t) = (r^*_{1}(t), r^*_{2}(t), \cdots, r^*_{n}(t))$. Let $\theta^*(t) = \frac{\sum\limits_{k = 1}^nkp(k)I^*_{k}(t)}{\langle k \rangle}$. And there exist adjoint varibles $\lambda_{1k}(t)$, $\lambda_{2k}(t)$, $\lambda_{3k}(t)$, $\lambda_{4k}(t)$ that satisfy

with transversality condition

In addition, the optimal control $r^*_{k}(t)$ is given by

Proof. According to the Pontryagin Maximum Principle [cankaowenxian] with the Hamiltonianfunction, the adjoint equations can be determined by the following equations

Furthermore, by the necessary condition, we have

So the optimal control problem can be determined by the following system:

with initial condition

and transversality condition

6.   Stochastic and numerical simulations

Owing to our models are network-based models, it is indispensable to carry out stochastic simulations as well as numerical simulations. The validity of the model can be better verified by means of combinating numerical simulations of differential system and node-based stochastic simulations.

6.1.   Experimental setup

In the following examples, numerical simulations is implementedby using the classic BA scale-free network model generation algorithm. We consider scale-free networks with degree distribution $P(k) = \eta k^{-2.5}$ for $h = 1, 2, ..., 100$, where the constant $\eta$ is chosen to maintain $\sum\limits_{k = 1}^{100}p(k) = 1$. Under the aforementioned conditions, we can obtain $\langle k \rangle = 1.7995$, $\langle k^2 \rangle = 13.8643$, $\langle k^3 \rangle = 500.7832$.

Optimality system $(5.9)$ is solved by using the forward-backward Runge-Kutta fourth order method ^[36].

he setting of the parameters of the models are summarized in Table 1. To note that the parameters can be changed according to various application scenarios.

6.2.   Experimental results

Example 1. The first example is used to verify the threshold and stability results of the model. For the network, in the case of Figure 3 (a), (b), (c), we draw the trajectory of the average number of infected $\sum\limits_{k = 1}^{n}p(k)I^*_{k}$ when $R_{0} < 1$, $R_{0} = 1$, and $R_{0} > 1$.

Example 2. In the second example, we study the case of forward bifurcation. We choose parameter $\alpha = 22, 25, 28, 31$ to describe the bifurcation diagram in Figure 3, showing a positive bifurcation when $R_{0} = 1$. An interesting observation is that the local balance curve has an "S-shaped" $\alpha$ value. In this case, there is backward bifurcation from an endemic equilibrium at a certain value of $R_{0}$, which leads to the existence of multiple characteristic balances.

Example 3. The third example shows the effectiveness of the optimal control strategy. We choose parameter $\alpha = 31$. The initial conditions are $S_{k}(0) = 0.9$, $E_{k}(0) = 0$, $I_{k}(0) = 0.1$, $R_{k}(0) = 0$. For the optimal control problems, the weights parameters are set as $A_{k} = 0.8$, $T = 10$. Through strategy simulation, the mean value of optimal control is obtained $\langle r \rangle = \frac{1}{T}\int_{0}^{T}r^*_k(t)dt$. The average value of the control measures is shown in the middle column of Figure 4. The left column shows the average number of infections for different control strategies, and the right column shows the cost. It can be seen from Figure 3 that in each case, the optimal control effect has reached the maximum control effect, but the cost is lower.

7.   Conclusions

In order to understand the impact of limited medical resources and population heterogeneity on disease transmission, this paper proposes a SEIR model with a saturated processing function based on a complex network. The model proved that there is a threshold $R_{0}$, which determines the stability of the disease-free equilibrium point. In addition, we also proved that there is a backward branch at $R_{0} = 1$. In this case, people cannot eradicate the disease unless the value of $R_{0}$ decreases such that $R_{0} < \hat{R}_{0} < 1$. In summary, in order to control or prevent infectious diseases, our research results show that if a backward bifurcation occurs at $R_{0} = 1$, then we need a stronger condition to eliminate the disease. However, if a positive bifurcation occurs at $R_{0} = 1$, the initial infectious invasion may need to be controlled at a low level, so that disease or death or close to a low endemic state of stability.

For the SEIR model with saturation processing function, in order to formulate effective countermeasures, two methods of constant control and time-varying control are introduced to control the model. For the optimal defense situation, the optimal control problem is discussed, including the existence and uniqueness of the optimal control and its solution. Finally, we performed some numerical simulations to illustrate the theoretical results. Optimal control does achieve a balance between control objectives and control costs.

Acknowledgments

Research Project Supported by National Nature Science Foundation of China (Grant No. 1191278), the Fund for Shanxi 1331KIRT and Shanxi Natural Science Foundation (Grant No. 201801D121153, 201901D111179).

Conflict of interest

The authors declare there is no conflict of interest in this paper.