Dynamic stability and optimal control of SIS q I q RS epidemic network.

We develop a complex network-based SIS q I q RS model, calculate the threshold R 0 of infectious disease transmission and analyze the stability of the model. In the model, three control measures including isolation and vaccination are considered, where the isolation is structured in isolation of susceptible nodes and the isolation of infected nodes. We regard these three kinds of controls as time-varying variables, and obtain the existence and the solution of the optimal control by using the optimal control theory. With regard to the stability of the model, sensitivity analysis of the parameters and optimal control, we carry out numerical simulations. From the simulation results, it is obvious that when the three kinds of controls exist simultaneously, the scale and cost of the disease are minimal. Finally, we fit the real data of COVID-19 to the numerical solution of the model.

Introduction

The spread of infectious diseases has been a major threat to human health for a long time, because large-scale population loss can occur when the epidemic breaks out. A large number of scholars have studied infectious diseases from various aspects, especially in infectious disease control. Several scholars have analyzed and studied network-based epidemic models [1], [2], [3], [4], [5], [6], ordinary differential models and fractional order derivative models [7], [8], [9], [10], [11], [12]. Research on traditional infectious disease models is relatively mature, but in these models, individuals are considered to be uniformly mixed, without considering the relationship between individuals. The number and probability of contact between each individual and other individuals are different, thus forming a network. Considering the spread of infectious diseases based on such a network will be closer to reality. Some scholars have studied scale-free networks or small-world networks by considering them into the infectious disease model, and many important results have been obtained. Li et al. [1] developed a nonlinear SIQS model on the network. They considered the problem of optimal isolation control, and studied the optimal control under different network structures by numerical simulation. Zhang et al. [3] proposed the SIQS model on the network, considered three control measures in the model, and discussed the time-varying optimal control problem of the model. Xu et al. [4] developed a SIVRS model based on complex networks considering virus variation factors and formulated an optimal control problem to reduce the spread of the virus. Liu et al. [6] studied the spread of epidemics on networks using the SIR pairwise approximation model, considering the optimal vaccination control problem, and simulated them on different networks.

There are numerous different types of infectious diseases, each of which is related and distinct, so many scholars have conducted research on a particular disease, such as COVID-19 [13], [14], [15], [16], Malaria [17], [18], and other infectious diseases [19], [20]. It is extremely important to control these infectious diseases. Due to the rapid spread of COVID-19, the epidemic situation is accelerating. On March 11th, 2020, the World Health Organization (WHO) officially declared that the infectious disease constitutes a global epidemic. COVID-19 outbreaks are fast and rapid, and there are many infectious disease models that have been developed and the magnitude of their outbreaks is worth studying. Many models have been developed based on disease transmission, some of which are related to estimating the risk of virus transmission and predicting the epidemic of disease [21], [22], [23], and some of which are related to the influence of media reports, control and travel restrictions on disease [24]. Wang et al. [16] developed an infectious disease model based on the spread of COVID-19 on a complex network. Isolation of susceptible individuals was included in the model, and data analysis was used to determine when 16 cities, including Wuhan, could resume work. Tchoumi et al. [18] developed a dual Malaria-COVID 19 model and considered the problem of optimal control of the model. Through two control measures, the number of people infected with these two diseases decreased. Tang et al. [21] developed a compartmental model based on the transmission of COVID-19 and interventions (including isolation of susceptible individuals), and calculated the basic reproduction number using two methods. From the numerical simulation, they concluded that it is effective to track and isolate close contacts.

The control of infectious diseases is extremely important. Many scholars have studied optimal control, including isolation, vaccination and treatment, aiming at reducing the scale of disease and controlling costs [25], [26], [27], [28], [29], [30]. Abboubakar et al. [25] developed a new model of vector-borne infectious disease and introduced four control functions in the model. In numerical simulations, they used the available data of Chikungunya to calibrate the parameters of the model by the least square method to verify the feasibility of the model. Richard et al. [26] modeled an infectious disease with an age structure. Considering the lack of treatment and vaccines, they focused on control measures with non-pharmacological interventions and identified a control scheme to minimize costs. Zhou et al. [28] mainly considered the reaction–diffusion SIR epidemic system. They studied two control strategies: vaccine and treatment, and verified the results by numerical simulations.

It can be clearly seen from the above literature that some models are built on complex networks, and some models are aimed at a specific infectious disease. They include vaccination or isolation in the models to consider optimal control. In time of infectious disease outbreaks, not only are those who are infected isolated, but for those who are susceptible. Isolation measures such as home isolation are required to prevent being infected. However, in the existing complex network models, isolation control of susceptible individuals is not included. We think it is necessary to take this factor into account in the models. Isolation of susceptible individuals, such as travel restrictions, home isolation and city closures, is the most effective and fastest control measure before finding a treatment and vaccine. In [31], Yang et al. pointed out that if the closure of the city was delayed by 5 days during the COVID-19 outbreak, the total transmission scale would increase nearly three times. It is necessary to isolate susceptible individuals in time. Thus, we build an SISIRS model on complex networks. In this model, we consider the isolation of susceptible and infected nodes and vaccination measures. We calculate the threshold to measure whether the outbreak occurred, and analyze the global asymptotic stability of the disease-free equilibrium and the endemic equilibrium. Besides proving the existence and the solution of optimal control, we also compare different control strategies to make the scale of the infectious disease as small as possible and the cost-effectiveness as good as possible. The numerical simulation consists of four parts. The stability of the model is verified, then the sensitivity of the parameters is analyzed, the simulation is performed from the perspective of optimal control, and finally, the numerical solution of the real data and the model is fitted. It can be observed that when three controls exist at the same time, the scale and cost of the disease are the smallest. Our network-based infectious disease model can be used for infectious diseases, such as SARS, A/H1N1 and COVID-19, which have high infection rates and spread rapidly, infect a large numbers of people in a short time, and in some cases lead to a global pandemic.

The following text consists of the following sections. In Section 2, we introduce the network-based SISIRS model. In Section 3, the expressions for the basic reproduction numbers are given and the stability of the model is analyzed. In Section 4, the optimal control of time-varying variable model is discussed. In Section 5, we validate the conclusions of the first two sections using numerical simulations and verify the validity of the model by fitting to the data. In Section 6, we summarize the full text.

Model formulation

Based on complex network, we establish an SISIRS model for epidemic disease. The status of each node is determined by the following five conditions: susceptible (S), infected (I), susceptible and isolated (S), infected and isolated (I), recovery (R). The spread of the disease is shown in Fig. 1. It is assumed that the probability of infection of susceptible nodes after each contact with infected nodes is , and the number of contacts is . When susceptible nodes are in closing contact with infected nodes but not infected, they are isolated at the rate of . Isolated susceptible nodes are desegregating at the rate and become susceptible nodes. After the infected nodes are diagnosed, it is required to be isolated at the rate of . The infected nodes and the isolated infected nodes become recovery nodes at the rates and , respectively. Recovery nodes become susceptible nodes by losing immunity at rate . The vaccination rate for susceptible nodes is .

Fig. 1

The state-transition rules of the SISIRS model.

According to the above-mentioned disease transmission process, our model is as follows: where , , , and represent the relative densities of susceptible nodes, infected nodes, isolated susceptible nodes, isolated infected nodes and recovery nodes of degree respectively, here . represents the probability that any link points to the infected node, where is the average degree and is the degree distribution.

The set defined by , , is a compact attracting positively invariant of model (1). Obviously, the disease-free equilibrium of model (1) is

Equilibria and global stability

Before studying the dynamic behavior of the model (1), we calculate the basic reproduction number of the model by the next generation matrix method [32], which is the key to measure whether the disease broke out or not.

The right-hand side of model (1) could be written as where

Let vector and vector . The Jacobian matrix of and at the disease-free equilibrium is given by where

Thus, the basic reproduction number can be calculated as , where represents the spectral radius of the matrix.

If , then the disease-free equilibrium of model (1) is globally asymptotically stable.

Due to , we can consider the following auxiliary system:

From is compact attracting positively invariant set and model (2), we can get:

Next, we consider:

The Jacobian matrix of the model (4) at the disease-free equilibrium is as follows: where , , .

Obviously, has eigenvalues that are , eigenvalues that are , and the remaining eigenvalues are the same as those of the matrix . According to [32], is locally asymptotically stable and , where represents the maximum value of real part of eigenvalue. The solutions of the model (4) tend to as . For model (1), it can be obtained by comparison principle that , , , , . It implies that is globally attractive. Thus, is globally asymptotically stable. □

If and only if , model (1) has a unique endemic equilibrium .

Let the right hand of model (1) be zero and according to , we can obtain where , The following self-consistent equation can be obtained by bringing into : where , , .

Obviously, and . From (5), we can obtain: If , then there is a nonzero positive solution for (5). That is, if , then the model (1) has a unique endemic equilibrium . □

If , the endemic equilibrium of model (1) is globally asymptotically stable.

We consider auxiliary system (2), and prove the theorem based on Theorem 3.1 and Corollary 3.2 of [33]. Let , , . Let the right hand of the model be , which is a continuously differentiable cooperative map. And the following conditions are satisfied: (i) , and is bounded for all with ; (ii) is irreducible for every ; (iii) choose any , for , for , for ,

for ,

Therefore, the model (2) is strictly sublinear on . We know that when , , is globally asymptotically stable by Corollary 3.2 of [33]. □

Optimal control

With the outbreak of COVID-19 epidemic in 2019, we know that it is necessary to isolate susceptible and infected individuals. In the early stages of an outbreak, there is no definite treatment method, and the most effective strategy to control the epidemic is isolation. Of course, another control strategy for susceptible individuals is vaccination. Therefore, in this section, we use the isolation rate for susceptible individuals, the isolation rate for infected individuals and the vaccination rate as time-varying control variables, and then discuss the optimal control problem of model (1) by using the optimal control theory.

Let be a time-varying control variable. And the control set is , where is the terminal time. The model (1) can be written in the following form:

In order to reduce the density of infected persons by control, we consider the following objective function: where Lagrangian , and represent different weights. Thus, the optimal control problem is to find a such that .

There exists an optimal control such that , subject to controlled model (6) of state variables.

We prove it using Cesari Theorem [34]. The following five conditions are satisfied: (i) it is easy to know the control space is a compact and closed; (ii) according to the continuum theory of differential systems [35], for any control variable , the solution of model (6) exists; (iii) let the right-hand of the model (6) be , and , . By is bounded by the linear combination of ; (iv) , , , thus the function is concave in ; (v) since , Let , such that . Therefore, according to Theorem 4.1. and Corollary 4.1. of [34], we can get the existence of the optimal solution of model (6). □

By Pontryagin’s principle of [34], define the Hamiltonian function: where are the adjoint variables. Next, we give the solution to the optimal control problem.

Let be the optimal solution of model (6) under optimal control , where . The adjoint variables satisfy: with transversality condition . The optimal control is as follows:

According to the Pontryagin’s minimum principle of [34], we can obtain the following dynamics of the adjoint variables: with transversality condition . After a simple calculation we can get (7). Next, under the optimality condition, we obtain Thus, the optimal control of model (6) whose compact expression is given by (8). □

Numerical simulations

We simulate the model (1) and the model (7) on a scale-free network (BA network) with power-law distribution, and take . The average densities of the five states are , , , and respectively, where and . All parameters are non-negative.

At first, we verify the conclusion of the third section by numerical simulation. Fig. 2 shows the dynamics of susceptible nodes, infected nodes, isolated susceptible nodes, isolated infected nodes and recovery nodes. Fig. 2(a) illustrates the global asymptotic stability of the disease-free equilibrium point, taking the following parameter values and initial values: , , , , , , , , and , where represents the initial value of the node with degree . Fig. 2(b) indicates the global asymptotic stability of the endemic equilibrium point. The values of the parameters and the initial values taken are as follows: , , , , , , , , and .

Fig. 2

The dynamics of the susceptible nodes, infected nodes, isolated susceptible nodes, isolated infected nodes and recovery nodes.

From the expression of , the relationship between each parameter and can be easily obtained. Next we perform sensitivity analysis on the parameters so as to observe how they affect the disease transmission. In Fig. 3, we simulate the effects of infection rate and cure rate , vaccination rate and isolation rate of infected nodes on , respectively. High infection rates are shown in Fig. 3 to be critical for disease transmission and to reduce the threshold for disease outbreaks when coupled with control measures.

Fig. 3

The relationship between with and , and , and respectively.

In Fig. 4, we simulate the effects of and , and immune loss rate , and on , respectively. In Fig. 4(a), it is shown that vaccination and isolation together significantly decrease the value of ; in Fig. 4(b), it is shown that even though vaccination reduces , it has a modest effect on the increase of with the loss of individual immunity. Compared with Fig. 4(b), the cure rate in Fig. 4(c) has a greater influence on the decrease of when immune loss occurs.

Fig. 4

The relationship between with and , and , and respectively.

Next, we use the iterative method of “forward–backward scanning” [36], [37] to numerically solve the optimal control model, in which the solver is the fourth-order Runge–Kutta algorithm. Here, , , , , , . Set the weights parameters as follows: , and . The initial conditions are , , , , . Define the average control strength as: , , . Then, we simulate the cases with one control variable, two control variables and three control variables respectively. We compare the effects of uncontrolled, constant and optimal control on the average infected densities. When there is one control variable or two control variables, the values of the other control parameters that are constants are shown in Table 1, where “-” indicates that the parameter is taken as a variable.

Table 1

The values of the parameters and .

Parameter	δ_k	q_S,k	q_I,k	Parameter	δ_k	q_S_k	q_I,k
Figs. 5(a) and 6(a)	–	0.02	0.1	Figs. 7(a) and 8(a)	–	–	0.1
Figs. 5(b) and 6(b)	0.28	–	0.1	Figs. 7(b) and 8(b)	–	0.02	–
Figs. 5(c) and 6(c)	0.28	0.02	–	Figs. 7(c) and 8(c)	0.28	–	–

When there is only one control variable, Fig. 5 shows effect of control variables as , and of the disease, respectively. In Fig. 5, the infection density is highest in the absence of control, and the disease size is smallest under optimal control. The value of the objective function shown in Fig. 6 indicates the control cost required for the control measures, corresponding to the three cases in Fig. 5. From Figs. 6(a) and 6(c), we can see that the control cost of optimal control is lower than that of the other cases, but the cost of optimal control is not the lowest in Fig. 6(b).

Fig. 5

The dynamics of the infected nodes when the control variables are , and , respectively.

Fig. 6

The values of objective function when the control variables are , and , respectively.

When two control variables are considered, there are three combinations. Similarly, we compare the effects of uncontrolled, constant and optimal control on the average infected densities. Fig. 7 tells us that among the three combinations, the infected densities under optimal control are the lowest. The values of the objective function show in Fig. 8 correspond to the three cases in Fig. 7. Although optimal control can better control the epidemic, the objective function value under optimal control is not the minimum as showing in Fig. 8(a).

Fig. 7

The dynamics of the infected nodes when the control variables are and , and , and , respectively.

Fig. 8

The values of objective function when the control variables are and , and , and , respectively.

Therefore, we consider the situation that three kinds of controls exist at the same time, and compare the effects of minimum, average, maximum and optimal controls on the disease by simulation. Fig. 9(a) shows that the densities of infected nodes under optimal control and maximum control is lower than that under minimum control and average control. However, as showing in Fig. 9(b), the control cost is the lowest under the optimal control. Fig. 9(c) tells us that with the passage of time, the average control intensity becomes smaller, which is the reason why the control cost is the lowest under optimal control.

Fig. 9

Control variables , and exist simultaneously. (a) The dynamics of the infected nodes. (b) The values of objective function . (c) Average strength of control variables.

We apply the SISIRS model to real examples to show the validity and accuracy of this model. 2019 COVID-19 is used as example for numerical simulation and fitting. The parameters in the model are estimated using the Markov Chain Monte Carlo (MCMC) method [38]. The sampling method is the Metropolis–Hastings (M–H) algorithm, and the convergence is evaluated using the Geweke convergence diagnostic method during the iterative process. The number of COVID-19 cases is taken from the official web site of the National Health Commission of the People’s Republic of China [39]. During the COVID-19 outbreak, we select the time period from the start of mass nucleic acid testing across China (February 12, 2021) to the end of the Wuhan closure (April 7, 2021) because the number of reported cases was inaccurate before mass nucleic acid testing was performed. The values and sources of the parameters and initial values are shown in Table 2.

Table 2

Values and sources of parameters and initial values.

Parameter/Initial value	β	c	q_S	α	q_I	γ_I	γ_q
Value	0.5	1.6897	0.51907	0.014733	0.14323	1/7	1/14
Source	[40]	MCMC	MCMC	MCMC	MCMC	[40]	[40]
Parameter/Initial value	σ	δ	S(0)	I(0)	S_q(0)	I_q(0)	R(0)
Value	0.0070164	0.43313	0.9651	0.0002	0.0344	0.0003	0
Source	MCMC	MCMC	[39]	MCMC	MCMC	[39]	[39]

Fig. 10 shows the fit of the actual number of cases and the numerical solution of the model for COVID-19. As seen in the figures, the degree of fit between the actual and simulated numbers can indicate the validity of our model.

Fig. 10

Number of cases of COVID-19 (February 12, 2021–April 7, 2021). The red circles indicate the actual number of cases and the solid black line indicates the numerical solution of the model.

Conclusion and discussion

Through years of observation on the spread of infectious diseases, we can find that for some infectious diseases, it is necessary to isolate susceptible populations, such as A/H1N1 and COVID-19. During the outbreak of the disease, especially during the COVID-19 outbreak, Wuhan, China took measures to close the city, while other regions also implemented strategies such as closing villages and communities to control people’s travel. People are asked to be isolated at home, and they should not go out unless necessary. Other measures to prevent transmission to susceptible individuals include vaccination.

Accordingly, in addition to vaccination and isolation of infected nodes, we consider isolation of susceptible nodes into the model, and built an SISIR epidemic model on the network. The analysis show that when , the disease-free equilibrium is globally asymptotically stable; When , there exists a unique globally asymptotically stable endemic equilibrium for the model (1). Next, we consider the three control as time-varying variables, and prove the existence and solution of optimal control by using the optimal control principle and Pontryagin’s principle. Simulations are carried out on the BA network, which verifies the previous conclusion. At the same time, we also simulated the situation that there is one control, two controls and three controls. Simulation results show that the density of infected nodes is lowest and the cost is lowest when all three controls are present simultaneously, resulting in optimal control. Finally, the numerical solution of the model is fitted with the real data of COVID-19, which shows that our model is effective. There are numerous kinds of complex networks. Using the model, we can also consider the influence of different topological structures on infectious disease control. Moreover, the control costs of non-pharmaceutical and pharmaceutical interventions can be compared to develop optimal strategies.

CRediT authorship contribution statement

Xinjie Fu: Writing – original draft, Writing – review & editing. JinRong Wang: Writing – original draft, Writing – review & editing, Supervision.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.