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

Abstract: In order to study the impact of limited medical resources and population heterogeneity on disease transmission, a SEIR model based on a complex network with saturation processing function is proposed. This paper first proved that a backward bifurcation occurs under certain conditions, which means that R0 < 1 is not enough to eradicate this disease from the population. However, if the direction is positive, we find that within a certain parameter range, there may be multiple equilibrium points near R0 = 1. Secondly, the influence of population heterogeneity on virus transmission is analyzed, and the optimal control theory is used to further study the time-varying control of the disease. Finally, numerical simulations verify the stability of the system and the effectiveness of the optimal control strategy.


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 .
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] h(I) = 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, α 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.

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]. θ(t) describes a link pointing to an infected node, which satisfies the relation kp(k) describes the average degree. Here, we assume that the connectivities of nodes in network at each time are uncorrelated. rI k (t) 1+αθ(t) represents the recovery of the k-th infection group after treatment.
1+αθ(t) = rθ(t) 1+αθ(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: (i) all new births are susceptible; (ii) the total number of nodes is constant, and the number of deaths is equal to the number of births, so π = µ; (iii) 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 k N k = N.
For the practice, the initial condition for model (2.1) satisfy: The probability 0 ≤ θ(t) ≤ 1 describes a link pointing to an infected host, which satisfies the relation At t time, the global ensity is 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

Equilibria and basic reproduction number
n} is a positively invariant for model (2.1).
Proof. Assuming that E = (S 1 , E 1 , I 1 , R 1 , · · · , 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 This implies that has a positive solution. While Then k . 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.

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 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×2n is given by: 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.
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 π = µ, we have Because 0 < I k (t) < 1 for all k, we know that 0 < θ(t) < 1. Therefore The equation holds if and only if I k (t) = 0. For the limit system E k (t) = −(µ + δ)E k (t), it is easy to see that lim t→+∞ E k (t) = 0, For the limit system R k (t) = −µR k (t), it is easy to see that lim t→+∞ R k (t) = 0.Finally, from S k (t)+ E k (t)+ I k (t)+R k (t) = 1, we get lim t→+∞ 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 whenR 0 < 1. In addition, you can see that R 0 <R 0 . This prompts us to think about whether the disease can persist even if R 0 < 1.

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 θ > 0. From F(θ) = 0, we can see that the unique equilibrium should satisfy the following equation : The denominator and numerator are multiplied by the quantity δ k 2 (µ+δ)(r+µ+γ+µ d ) k , which can be expressed by R 0 and θ as follows: If the local θ is a function of R 0 , the sign of the bifurcation direction is the slope at (R 0 , θ) = (1, 0) (cf. Figure 2). More specifically, if the derivative is positive at the critical value (R 0 , θ) = (1, 0), that is ∂θ ∂R 0 (R 0 ,θ)=(1,0) > 0 the endemic equilibrium curve diverges forward. Conversely, if the derivative is negative at the critical value (R 0 , θ) = (1, 0), that is ∂θ ∂R 0 (R 0 ,θ)=(1,0) < 0 We can from Eq (4.2) that To conclude, we have the following theorem: It can be seen from Theorem 4 that the nonlinear processing function does play a key role in causing backward bifurcation. More precisely, if α is large enough to satisfy condition (4.3), backward bifurcation will occur. Otherwise, when this effect is weak, there is no backward bifurcation.

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 Our objective is to find a optimal control r * k (t) such that T ]} is the control set. Next we will analyze the optimal control problem.

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 ∈ Usuch that the constraint dynamical system is solvable. (C3) F(X(t), r k (t)) is bounded by a linear function in X.
Proof. Let us show that the five conditions in Lemma 3 hold true.
i) 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 , we have r ∈ (L 2 [0, T ]) n . The closedness of U follows from 0 ≤ r = lim n→∞ r n ≤ 1.
Let r 1 , r 2 ∈ U, η ∈ (0, 1), we have 0 ≤ (1 − η)r 1 + ηr 2 ≤ 1 ∈ (L 2 [0, T ]) n as (L 2 [0, T ]) n is a real vector space. ii) For any control variable r ∈ U, the solution of system (12) obviously exists following from the Continuation Theorm for Differential Systems [34]. iii) 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.

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 λ 1k (t), λ 2k (t), λ 3k (t), λ 4k (t) are the adjoint variables to be determined later.

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. 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.

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   Example 2. In the second example, we study the case of forward bifurcation. We choose parameter α = 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" α 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 α = 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 r = 1 T T 0 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.

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 <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.