Graph Propagation Dynamics of an Epidemic Model with Heterogeneous Control Strategies on Complex Networks

1. Introduction

The prevalence and outbreak of infectious diseases is always an important issue affecting the national economy and people’s livelihood [1]. Researches on the transmission mechanism and dynamics of infectious diseases have attracted the attention of many scholars. Specifically, mathematical models have been a potent tool for forecasting the trajectory of infectious diseases and evaluating different prevention and control measures [2,3]. Through the application of optimal control theory and qualitative assessment of mathematical models, the effective management, reduction, and potential eradication of an infectious disease can be achieved. [4,5,6,7,8]. Based on different mechanisms of disease transmission, Kermack and McKendrick introduced the SIS and SIR epidemic models [9,10,11]. Subsequent to their work, a variety of epidemic models have emerged to shed light on the progression of various diseases and to offer insightful control measures.

Any research progress on the dynamics and control of epidemic transmission may have a significant impact on the prevention and control of infectious diseases, which has attracted widespread attention from scientific research [12,13]. During the transmission of infectious diseases, the infected transmission mainly occurs through the contact between susceptible and infected individuals, then quarantine has become a commonly employed strategy for controlling disease. For example, Hethcote et al. studied SIQS and SIQR epidemiology models with different incidence [14]. The Refs. [15,16,17] have also conducted qualitative analysis of mathematical models with infected isolation on the assumption of uniformly mixing population. Compared with the classics propagation models, heterogeneous networks own a better reflection with the statistical properties of real social networks. Li et al. introduced an SIQRS epidemic model based on scale-free networks, and their theoretical results indicate that the epidemic threshold is closely linked to the network’s structural characteristics [18]. Kang et al. explored how the quarantine based on node degree, as well as constant rates of birth and death, affect disease dynamics on scale-free networks by utilizing a novel deterministic model [19]. Chen et al. formulated SVIQR and SVIQS models with heterogeneous quarantine rate and immunity rate to investigate the effects of intervention methods on the behavior of infectious diseases within heterogeneous networks [12]. Cheng et al. considered a network-based SIQS model with non-monotone incidence rate and studied the global stability of equilibria [20]. Zhao et al. developed a stochastic SIQS epidemic model on scale-free networks, taking into account random fluctuations in infection rates [21]. However, those models either ignored the demographics or assumed that the recruitment of population is constant.

Consistently implementing control measures can effectively reduce the scale of infection, however, if considering the economic factors, relying solely on continuous control measures may not be sufficient to achieve optimal control goals when the disease prevalence is substantial. Therefore, A viable approach is to employ time-varying control within certain bounds to strike a balance between control objectives and the associated costs via optimal control theory. Optimal control applications in epidemic dynamics primarily concentrated on homogeneous contact networks [22,23,24,25,26,27]. Recently, many intriguing works on the optimal control of epidemics on heterogeneous networks have emerged [18,28,29,30,31,32,33,34,35]. Li et al. presented a nonlinear SIQS epidemic model on networks and explored the issue of optimal quarantine control to minimize the cost of control measures [32]. Zhang et al. examined the optimal control of an SIQRS model that includes vaccination on a network, and studied the effects of different control strategies [34]. Yang et al. investigated the stability and optimal control on treatment for SIS epidemic systems combining in directed networks [35]. However, the above results only considered the optimal control strategy for isolation and did not compare it with other forms of heterogeneous isolation strategies.

The above works provide strong theoretical supports for the prevention of infectious diseases through isolation measures. However, appropriate vaccination for susceptible individuals is also crucial for preventing and controlling infection. On the basis of above consideration, we study an SIS epidemic model with heterogeneous immune for susceptible individuals and quarantine for infected one, which also consider the demographics of individuals on heterogeneous networks. The organization of this paper is as follows. In Sect.2, we introduce an SIQS model with heterogeneous control strategies1 of individuals on complex networks and provide some description and assumptions. In Sect.3, we obtain the equilibria and basic reproduction number. The main results are showed in Sect.4, we establish uniform persistence and analyze the global stability of the disease-free and endemic equilibria of system by different mathematical methods. The optimal time-varying quarantine control for the model is considered in Sect.5. In Sect.6, some numerical simulations are conducted to illustrate and supplement the analysis results. Finally, a brief conclusion is given in Sect.7.

2. Description and formation of epidemic models

Consider a population situated on a complex network N, in which each node of N represents either an individual or an unoccupied state. Every node has three optional states: susceptible state(S), quarantined state(Q), infected state(I) [21]. We divided the population into n groups, let $S_k\left(t\right)$, $I_k\left(t\right)$, $Q_k\left(t\right)$ denote the densities of susceptible, infected and quarantine nodes (individuals) with connectivity (degree) k at time t, respectively. The mechanism of our model mainly involves three factors. The schematic diagram of the model is shown in Figure 1, the meaning of parameters is listed in Table 1.

(1) Birth and Death: Each vacant node i in the network randomly selects a neighbor. If the neighbor is a vacant node, the state of i remains unchanged. If the neighbor is a non-vacant node, the vacant node i will be activated to generate a new susceptible nodes at the birth rate b. All non-vacant nodes become vacant nodes at a natural death rate d per unit time. We assume that each non-vacant node has the same birth contact ability A (where $A=1$) due to physiological constraints.

(2) Immunization and Quarantine ($I\to Q$): At each time step, susceptible individuals with degree k are immunized with the immune rate ${\delta }_k$. The infected nodes with degree k will be quarantined with rate ${\beta }_k$. The quarantined individuals will recover to susceptible node with rate $\eta$. Nodes with the same degree have identical quarantine (immunization) strategies, while those with different degrees have different strategies.

(3) SIS epidemic framework: Infection $S\to I$: At the initial moment, randomly select some nodes as infected nodes. At each time step, the possibility that each infected node i connect to its neighboring nodes is $\frac{\phi \left(k_i\right)}{k_i}$, where $\phi \left(k\right)$ represents the infectivity of infected nodes with degree k, and $\phi \left(k\right)$$=k$ [36,37,38], $\phi \left(k\right)$=A [39], $\phi \left(k\right)$=$k^m$ [40], $\phi \left(k\right)$=$\frac{a{k}^m}{1+b{k}^m}$ [41]. If an infected node i interacts with a susceptible node j along a connecting edge, node j has a possibility of being infected by i at a transmission rate ${\lambda }_{ij}$. For the node with degree k, the overall transmission rate is $\lambda \left(k\right)$.

Recovery $I\to S$: Each infected node reverts to being susceptible at recovery rate $\gamma$.

According to above assumption, the model is as follows,

$$\left\{\begin{array}{c}\frac{d{S}_k\left(t\right)}{dt}=bk[1-N_k\left(t\right)]{\Psi }_k-d{S}_k\left(t\right)-\lambda \left(k\right)(1-{\delta }_k)S_k\left(t\right)\Theta \left(t\right)+\gamma I_k\left(t\right)+\eta Q_k\left(t\right),\hfill \\ \frac{d{I}_k\left(t\right)}{dt}=\lambda \left(k\right)(1-{\delta }_k)S_k\left(t\right)\Theta \left(t\right)-({\beta }_k+\gamma +d)I_k\left(t\right),\hfill \\ \frac{d{Q}_k\left(t\right)}{dt}={\beta }_k{I}_k\left(t\right)-(\eta +d)Q_k\left(t\right).\hfill \end{array}\right.$$

(1)

with initial conditions $0\le S_k\left(0\right),I_k\left(0\right),Q_k\left(0\right)\le 1,k=1,2,\dots ,n$, where $\Theta \left(t\right)={\sum }_{i=1}^{n}p\left(i\right|k)\frac{\phi \left(i\right)}{i}{I}_i\left(t\right)$ stands for the probability that the infected contact through a link which originates from a node with degree k and points to the infected node. The ${\Psi }_k\left(t\right)={\sum }_{i=1}^{n}p\left(i\right|k)\frac{A}{i}{N}_i\left(t\right)$ denotes the probability that a node of degree k has a reproductive contact with an individual neighbor, the $\frac{1}{i}$ denotes the probability that an individual neighbor of vacant node with degree i. On uncorrelated networks, the likelihood is deemed to be independent of the connectedness of the originating node of the link. Thus, we note $\Theta \left(t\right)=\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)I_i\left(t\right)$, ${\Psi }_k\left(t\right)=\frac{1}{\langle k\rangle }{\sum }_{i=1}^{n}p\left(i\right)N_i\left(t\right)$. The average density of susceptible, infected and quarantined nodes is $S\left(t\right)={\sum }_{k=1}^n{S}_k\left(t\right)p\left(k\right)$, $I\left(t\right)={\sum }_{k=1}^n{I}_k\left(t\right)p\left(k\right)$, $Q\left(t\right)={\sum }_{k=1}^n{Q}_k\left(t\right)p\left(k\right)$.

Let $N_k\left(t\right)$ denotes the density of individuals of degree k at time t, and $N_k=S_k+I_k+Q_k.$ We develop the equation for $N_k$ as,

$$\frac{d{N}_k\left(t\right)}{dt}=b(1-N_k){\Psi }_k\left(t\right)-d{N}_k\left(t\right).$$

(2)

According to Ref. [42], we get ${lim}_{t\to \infty }{N}_k\left(t\right)=0$ when $b\le d$, no other dynamic patterns are present; when $b>d$, ${lim}_{t\to +\infty }{N}_k\left(t\right)=N_k{\left(t\right)}^{*}$, where $N_k^{*}=\frac{bk{\Psi }^{*}}{d+bk{\Psi }^{*}}$, ${\Psi }^{*}=\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\frac{ip\left(i\right)b{\Psi }^{*}}{d+bi{\Psi }^{*}}$.

Given that the original system and the limiting system exhibit identical long-term dynamic characteristics. In order to better analyze the stability of the model, we can only consider the limit system of (1) under the condition of $b>d$ and $N_k^{*}=S_k\left(t\right)+I_k\left(t\right)+Q_k\left(t\right)$ based on the analysis of $N_k\left(t\right)$. The corresponding limiting systems of model (1) are formulated as follows:

$$\left\{\begin{array}{c}\frac{d{S}_k\left(t\right)}{dt}=bk[1-N_k^{*}]{\Psi }^{*}-d{S}_k\left(t\right)-\lambda \left(k\right)(1-{\delta }_k)S_k\left(t\right)\Theta \left(t\right)+\gamma I_k\left(t\right)+\eta Q_k\left(t\right),\hfill \\ \frac{d{I}_k\left(t\right)}{dt}=\lambda \left(k\right)(1-{\delta }_k)S_k\left(t\right)\Theta \left(t\right)-({\beta }_k+\gamma +d)I_k\left(t\right),\hfill \\ \frac{d{Q}_k\left(t\right)}{dt}={\beta }_k{I}_k\left(t\right)-(\eta +d)Q_k\left(t\right).\hfill \end{array}\right.$$

(3)

with initial conditions $0\le S_k\left(0\right),I_k\left(0\right),Q_k\left(0\right)\le N_k^{*}<1,k=1,2,\dots ,n$.

The following proposition shows that the state space of the solutions of system (3) is positively invariant.

Proposition 1.

If ${(S_k\left(t\right),I_k\left(t\right),Q_k\left(t\right))}_k$ represents a solution to the system (3) that meets the initial conditions, then $0\le S_k\left(t\right)+I_k\left(t\right)+Q_k\left(t\right)\le N_k^{*}<1$, $k\in {\mathbb{N}}_n$ for any $t\ge 0$, i.e., Γ is positively invariant.

Proof. 

∀$x_0\in \Gamma$. Thus, $\forall t>0$, we have $N_{k0}\ge 0,\Psi \left(0\right)\ge 0$.

$$\frac{d\Psi \left(t\right)}{dt}=(b-d)\Psi \left(t\right)-\frac{b\Psi \left(t\right)}{\langle k\rangle }\sum _{i}ip\left(i\right)N_i\left(t\right).$$

We obtain that $\Psi \left(t\right)=\Psi \left(0\right)e^{b-d-\frac{b}{\langle k\rangle }{\Sigma }_{i}ip\left(i\right)N_i\left(t\right)}\ge 0$. By $N_k\left(t\right)=S_k\left(t\right)+I_k\left(t\right)+Q_k\left(t\right)$ and system (2), the evolution of $N_k\left(t\right)$ is derived by the following equation:

$$\frac{d{N}_k\left(t\right)}{dt}=bk[1-N_k\left(t\right)]\Psi -d{N}_k\left(t\right)\ge -[bk\Psi +d]N_k\left(t\right).$$

We have $N_k\left(t\right)\ge N_k\left(0\right)e^{-(bk\Psi +d)t}\ge 0$. Then, from system (2), we also acquire

$$\frac{d(1-N_k\left(t\right))}{dt}=-(bk\Psi +d)[1-N_k\left(t\right)]+d.$$

Obviously, $N_k\left(t\right)<1$, and $0\le S_k\left(t\right)+I_k\left(t\right)+Q_k\left(t\right)=N_k\left(t\right)<1.$

Then, we verify $I_k\ge 0$ for any $k=1,2,\dots ,n$. If not, because $I_k\left(0\right)\ge 0$, there is exist a $j\in \{1,2,\dots ,n\}$ and $t_0\ge 0$ satisfy

$$t_0=inf\{t|I_j\left(t\right)=0,\frac{d{I}_j\left(t\right)}{dt}<0\}.$$

It is evident that $I_j\left(t_0\right)=0,I_j^{\prime }\left(t_0\right)<0$ for any $t\in [0,t_0]$, $I_j\left(t\right)\ge 0$, we have

$$\frac{d{I}_j\left(t_0\right)}{dt}=\lambda \left(j\right)(1-{\delta }_j)S_j\left(t_0\right)\Theta \left(t_0\right)-({\beta }_j+\gamma +d)I_j\left(t_0\right)=\lambda \left(j\right)(1-{\delta }_j)S_j\left(t_0\right)\Theta \left(t_0\right)<0.$$

Because of $\Theta \left(t_0\right)=\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)I_i\left(t_0\right)\ge 0$, we obtain $S_j\left(t_0\right)<0$. On the other hand, for any $t\in [0,t_0]$, $Q_j\left(0\right)\ge 0$, we attain

$$\frac{d{Q}_j\left(t\right)}{dt}={\beta }_j{I}_j\left(t\right)-(\eta +d)Q_j\left(t\right)\ge -(\eta +d)Q_j\left(t\right).$$

That is to say $Q_j\left(t\right)\ge Q_j\left(0\right)e^{-(\eta +d)t}\ge 0$. When $t\in [0,t_0]$,

$$\begin{array}{ccc}\hfill \frac{d{S}_j\left(t\right)}{dt}& =& bj[1-N_j\left(t\right)]{\Psi }_j-d{S}_j\left(t\right)-\lambda \left(j\right)(1-{\delta }_j)S_j\left(t\right)\Theta \left(t\right)+\gamma I_j\left(t\right)+\eta Q_j\left(t\right)\hfill \\ & \ge & -d{S}_j\left(t\right)-\lambda \left(j\right)(1-{\delta }_j)S_j\left(t\right)\Theta \left(t\right).\hfill \end{array}$$

Then, we acquire $S_j\left(t_0\right)\ge S_j\left(0\right)e^{-{\int }_0^{t_0}d+{\lambda }_j(1-{\delta }_j)\Theta \left(t\right)\mathrm{d}\tau }\ge 0$. It is contradictory, hence, $I_k\left(t\right)\ge 0,$ for any k and $t\ge 0$. We also can get $S_k\left(t\right),Q_k\left(t\right)\ge 0$, and $S_k\left(t\right),I_k\left(t\right),Q_k\left(t\right)<1$. That completes the proof. □

3. Equilibria and Basic Reproduction Number

In this section, the basic reproduction number $R_0$ is determined, which is the average number of new infections generated by a single newly infected individual throughout the entire infection period [20]. Considering the equivalent system,

$$\left\{\begin{array}{c}\frac{d{I}_k\left(t\right)}{dt}=\lambda \left(k\right)(1-{\delta }_k)(N_k^{*}-I_k\left(t\right)-Q_k\left(t\right))\Theta \left(t\right)-({\beta }_k+\gamma +d)I_k\left(t\right),\hfill \\ \frac{d{Q}_k\left(t\right)}{dt}={\beta }_k{I}_k\left(t\right)-(\eta +d)Q_k\left(t\right).\hfill \end{array}\right.$$

(4)

Denote X is the space of functions, $x=(I_1\left(t\right),Q_1\left(t\right),\dots ,I_n\left(t\right),Q_n\left(t\right))\in X$ is the state space for system (4).

$$\Gamma =\{x\in X:0\le I_k\left(t\right),Q_k\left(t\right)<1,k=1,2,\dots ,n\}.$$

(5)

The Proposition 1 can show that $\Gamma$ is the positive invariant for system (4). All viable steady states of system (4) are satisfied the following equation,

$$\left\{\begin{array}{c}\lambda \left(k\right)(1-{\delta }_k)(N_k\left(t\right)-I_k\left(t\right)-Q_k\left(t\right))\Theta \left(t\right)-({\beta }_k+\gamma +d)I_k\left(t\right)=0,\hfill \\ {\beta }_k{I}_k\left(t\right)-(\eta +d)Q_k\left(t\right)=0.\hfill \end{array}\right.$$

(6)

In particular, there exists disease-free equilibrium of (4) $E^0=(0,0,\dots ,0))$, and the positive equilibrium $E^{*}={(I_k^{*},Q_k^{*})}_k$, where $E^{*}$ satisfies that,

$$Q_k^{*}=\frac{{\beta }_k}{\eta }{I}_k^{*},{\Theta }^{*}=\frac{1}{k}\sum _{i=1}^n\phi \left(i\right)p\left(i\right)I_i^{*},I_k^{*}=\frac{\lambda \left(k\right)(1-{\delta }_k)N_k^{*}{\Theta }^{*}}{\lambda \left(k\right)(1-{\delta }_k)(1+\frac{{\beta }_k}{\eta }){\Theta }^{*}+{\beta }_k+\gamma +d}.$$

The basic reproduction number $R_0$ can be determined according to the following theorem.

Theorem 1.

Define the basic reproduction number

$$R_0=\frac{1}{\langle k\rangle }\sum _{i=1}^n\frac{\phi \left(i\right)p\left(i\right)\lambda \left(i\right)(1-{\delta }_i)N_i^{*}}{{\beta }_i+\gamma +d}.$$

There always exists a disease-free equilibrium $E^0$. If and only if $R_0>1$, the system (4) has unique positive equilibrium point $E^{*}$ which satisfy that

$$Q_k^{*}=\frac{{\beta }_k}{\eta }{I}_k^{*},{\Theta }^{*}=\frac{1}{\langle k\rangle }\sum _{i=1}^n\phi \left(i\right)p\left(i\right)I_i^{*},I_k^{*}=\frac{\lambda \left(k\right)(1-{\delta }_k)N_k^{*}{\Theta }^{*}}{\lambda \left(k\right)(1-{\delta }_k)(1+\frac{{\beta }_k}{\eta }){\Theta }^{*}+{\beta }_k+\gamma +d}.$$

(7)

Proof. 

Eq. (6) implies that

$$I_k=\frac{\lambda \left(k\right)(1-{\delta }_k)N_k\Theta }{\lambda \left(k\right)(1-{\delta }_k)(1+\frac{{\beta }_k}{\eta })\Theta +{\beta }_k+\gamma +d},$$

(8)

Substituting (8) into $\Theta$. Then we get,

$$\Theta =\frac{1}{\langle k\rangle }\sum _{i=1}^n\frac{\phi \left(i\right)p\left(i\right)\lambda \left(i\right)(1-{\delta }_i)N_i\Theta }{\lambda \left(i\right)(1-{\delta }_i)(1+\frac{{\beta }_i}{\eta })\Theta +{\beta }_i+\gamma +d}\triangleq \Theta f(\Theta ),$$

(9)

where

$$f(\Theta )=\frac{1}{\langle k\rangle }\sum _{i=1}^n\frac{\phi \left(i\right)p\left(i\right)\lambda \left(i\right)(1-{\delta }_i)N_i}{\lambda \left(i\right)(1-{\delta }_i)(1+\frac{{\beta }_i}{\eta })\Theta +{\beta }_i+\gamma +d}.$$

(10)

Because of $f^{\prime }(\Theta )<0$ and $\underset{\Theta \to +\infty }{lim}f(\Theta )=0$, it is clear that Eq. (9) has the positive solution if and only if $f\left(0\right)=\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\frac{\phi \left(i\right)p\left(i\right)\lambda \left(i\right)(1-{\delta }_i)N_i}{{\beta }_i+\gamma +d}>1$. That is to say, when $R_0>1$, Eq. (6) has unique positive solution ${\Theta }^{*}$, namely, a unique positive equilibrium point $E^{*}$ of system (4) exists described by equation (7) if and only if $R_0>1$. The proof is completed. □

4. Stability Analysis for SIQS Model

4.1. Stability Analysis of Disease-Free Equilibrium

Firstly, it is proved that the disease-free equilibrium (DFE) $E_0$ is locally asymptotically stable by the same method in [43] by the following theorem.

Theorem 2.

If $R_0<1$, then the DFE of system (4) is locally asymptotically stable and unstable if $R_0>1$.

Proof. 

The Jacobian matrix of system (4) at $E_0$ is:

$$A=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right),$$

where $A_{12}$ is the zero matrix with n order, $A_{21}={\beta }_k{E}_n$, $A_{22}=-(\eta +d)E_n$, $E_n$ is identity matrix, and

$$A_{11}=\left(\begin{array}{cccc}{\xi }_1{\phi }_1{p}_1-(\gamma +{\beta }_1)& {\xi }_1{\phi }_2{p}_2& \dots & {\xi }_1{\phi }_n{p}_n\\ {\xi }_2{\phi }_1{p}_1& {\xi }_2{\phi }_2{p}_2-(\gamma +{\beta }_2)& \dots & {\xi }_2{\phi }_n{p}_n\\ ⋮& ⋮& \dots & ⋮\\ {\xi }_n{\phi }_1{p}_1& {\xi }_n{\phi }_2{p}_2& \dots & {\xi }_n{\phi }_n{p}_n-(\gamma +{\beta }_n)\end{array}\right),$$

where ${\xi }_i=\frac{\lambda \left(i\right)(1-{\delta }_i)N_i^{*}}{\langle k\rangle },$${\phi }_i{p}_i=\phi \left(i\right)p\left(i\right)$. Let $y=\{y_1,y_2,\dots ,y_n,y_{n+1},\dots ,y_{2n}\}\in \Gamma$.

Subsequently, the system (4) can be reformulated in a vector format as:

$$\frac{dy}{dt}=Ay+N\left(y\right).$$

(11)

Let $N\left(y\right)$ represent a column vector and $N\left(y\right)\le 0$,

$$N\left(y\right)=-\left(\begin{array}{c}\lambda \left(1\right)(1-{\delta }_1)(I_1+Q_1)\Theta \\ \lambda \left(2\right)(1-{\delta }_2)(I_2+Q_2)\Theta \\ ⋮\\ \lambda \left(n\right)(1-{\delta }_n)(I_n+Q_n)\Theta \\ 0\\ 0\\ ⋮\\ 0\end{array}\right).$$

We consider the eigenvalues of matrix A are determined by $|\rho E_n-A_{11}\left|\right|\rho E_n-A_{22}|=0$. It is obvious that the characteristic equation has n multiple roots $-(\eta +d)$ from $|\rho E_n-A_{22}|=0$. Other eigenvalues of matrix A is determined by $|\rho E_n-A_{11}|=0$. Through some elementary transformations, $|\rho E_n-A_{11}|=0$ changes to

$$\left|\begin{array}{cccc}\rho -{\xi }_1{\phi }_1{p}_1+({\beta }_1+\gamma +d)& -{\xi }_1{\phi }_2{p}_2& \dots & -{\xi }_1{\phi }_n{p}_n\\ \frac{{\xi }_1}{{\xi }_2}(\rho +{\beta }_1+\gamma +d)& \rho +({\beta }_2+\gamma +d)& \dots & 0\\ ⋮& ⋮& \dots & ⋮\\ \frac{{\xi }_1}{{\xi }_n}(\rho +{\beta }_1+\gamma +d)& 0& \dots & \rho +({\beta }_n+\gamma +d)\end{array}\right|=0.$$

Case 1 If $\rho \ne -({\beta }_i+\gamma +d),$ then we get

$$(\rho +{\beta }_1+\gamma +d)\left[1-\frac{1}{\langle k\rangle }\sum _{i=1}^n\frac{\lambda \left(i\right)(1-{\delta }_i)N_i^{*}\phi \left(i\right)p\left(i\right)}{\rho +{\beta }_i+\gamma +d}\right]=0.$$

(12)

Note $G\left(\rho \right)=\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\frac{\lambda \left(i\right)(1-{\delta }_i)N_i^{*}\phi \left(i\right)p\left(i\right)}{\rho +{\beta }_i+\gamma +d}$. Obviously, the solutions of equation (12) are determined by $G\left(\rho \right)=1$. Since $G^{\prime }\left(\rho \right)<0$, $G\left(0\right)=R_0$, $\underset{\rho \to +\infty }{lim}G\left(\rho \right)=0,$${lim}_{\rho \to -({\beta }_i+\gamma +d)}G\left(\rho \right)=\infty$, hence, when $R_0<1$, the solutions of equation $G\left(\rho \right)=1$ are negative. Therefore, the equation (12) only has negative solutions. Case 2 If ${\prod }_{i=2}^n(\rho +{\beta }_i+\gamma +d)=0$, since the sum of all eigenvalues is equal to the trace of the matrix. Therefore, when $R_0<1$, ${\rho }_1={\xi }_1\phi \left(i\right)p\left(i\right)-({\beta }_1+\gamma +d)<0$, and ${\rho }_i=-({\beta }_i+\gamma +d)<0$.

Furthermore, because $N\left(y\right)\le 0$, thus $\frac{dy}{dt}\le Ay$. Since all the eigenvalues of matrix A are negative, the equilibrium point $E^0$ of system (4) is locally asymptotically stable when $R_0<1$. □

Then, we get the global asymptotic stability of $E^0$.

Theorem 3.

When $R_0<1$, the disease-free equilibrium $E^0$ of system (4) is global asymptotic stability.

Proof. 

From system (3), we get

$$\begin{array}{ccc}\hfill \frac{d\Theta \left(t\right)}{dt}& =& \frac{1}{\langle k\rangle }\sum _{k=1}^n\phi \left(k\right)p\left(k\right)\lambda \left(k\right)(1-{\delta }_k)(N_k^{*}-I_k\left(t\right)-Q_k\left(t\right))\Theta \left(t\right)-({\beta }_k+\gamma +d)\Theta \left(t\right)\hfill \\ & \le & \frac{1}{\langle k\rangle }\sum _{k=1}^n\phi \left(k\right)p\left(k\right)\lambda \left(k\right)(1-{\delta }_k)N_k^{*}\Theta \left(t\right)-({\beta }_k+\gamma +d)\Theta \left(t\right)\hfill \\ & =& ({\beta }_k+\gamma +d)\left(\frac{1}{\langle k\rangle }\sum _{k=1}^n\frac{\phi \left(k\right)p\left(k\right)\lambda \left(k\right)(1-{\delta }_k)N_k^{*}}{{\beta }_k+\gamma +d}-1\right)\Theta \left(t\right)\hfill \\ & =& ({\beta }_k+\gamma +d)(R_0-1)\Theta \left(t\right).\hfill \end{array}$$

We obtain $\Theta \left(t\right)\le \Theta \left(0\right)e^{({\beta }_k+\gamma +d)(R_0-1)t}$, while $\Theta \left(t\right)\ge 0$ from the Proposition 1. Therefore, if $R_0<1$, we can acquire $\underset{t\to +\infty }{lim}\Theta \left(t\right)=0$, and $\underset{t\to +\infty }{lim}{I}_k\left(t\right)=0$. From the second equation of system (4), we get $\underset{t\to +\infty }{lim}{Q}_k\left(t\right)=0$. This demonstrates that $E^0$ of system (4) is globally stable when $R_0<1$. The proof is finished. □

4.2. Global Stability of Endemic Equilibrium

In this section, we seek to study the conditions for the uniform persistence of the disease which is important in proving global stability of endemic equilibrium $E^{*}$. The result of uniform persistence is showed in the following Lemma.

Lemma 1.

When $R_0>1$, the system (4) is permanent, i.e., there exists $\epsilon >0$, such that

$$\underset{t\to +\infty }{lim}inf{\{I_k\left(t\right),Q_k\left(t\right)\}}_{k=1}^n\ge \epsilon ,$$

where $(I_k\left(t\right),Q_k\left(t\right))$ is any solution of system (4) satisfying (5).

Proof. 

Based on the Theorem 4.2 of [42], we need to verify that all hypotheses for system (4) are satisfied. By setting $C=\Gamma$, Condition (1) is fulfilled. Condition (2) is evidently met. For condition (3), notice that $A_{11}^T$ is irreducible and $a_{ij}\ge 0$, whenever $i\ne j$, hence there exists an positive eigenvector $\tilde{h}=(h_1,h_2,...,h_n)$ of $A_{11}^T$ and its corresponding eigenvalue is $s\left(A_{11}\right)$($s\left(A\right)=R{e}_i{\lambda }_i$), when $R_0>1$, $s\left(A_{11}\right)>0$. Let $h_{n+1}=h_{n+2}=\dots =h_{2n}=0$, then, we note $h=(h_1,h_2,\dots ,h_{2n})$, and it satisfies $A^{T}h=\lambda h$. Let $r={min}_{1\le i\le n}\left(h_i\right)$, for any $y\in \Gamma$, one has $h·y\ge r{\displaystyle \sum _{i=1}^n}{y}_i\ge r\sqrt{{\displaystyle \sum _{i=1}^n}{y}_i^2}$. Thus, condition (3) is verified. Since every part of $N\left(y\right)$ is negative and $h\ge 0$, hence, condition (4) is satisfied. To confirm condition (5), we set $G=\{y\in \Gamma |h·N\left(y\right)=0\}$. If $y\in G$, then

$$\sum _{i=1}^n{h}_i(y_i+y_{n+i})\lambda \left(i\right)(1-{\delta }_i)\Theta =0.$$

Given that every term in above sum is non-negative, one has $y_i=0$ for all $i=1,2,\dots n$. Therefore, the only invariant set corresponding to (11) that is contained within G is $y=0$. Condition (5) is satisfied. Hereby, all conditions are satisfied, the proof of uniform persistence is completed. □

From Lemma 1, it is clear that if $R_0>1$, the infection will always exists. The following theorem is showed that the global stability of endemic equilibrium.

Theorem 4.

Suppose that $(I_k\left(t\right),Q_k\left(t\right))$ is a solution of (4) satisfying initial condition. If $R_0>1$, then $\underset{t\to +\infty }{lim}{I}_k\left(t\right)=I_k^{*}$, $\underset{t\to +\infty }{lim}{Q}_k\left(t\right)=Q_k^{*}$, that is to say, the endemic equilibrium $E^{*}$ is global attractivity.

Proof. 

Firstly, from Lemma 1 and $k\in \{1,2,\dots ,n\}$, we get there exist $0<\epsilon \ll 1,$ such that, as t is sufficiently large, $I_k\ge \epsilon$ and $Q_k\ge \epsilon$ is satisfied. Hence, when $t>0$, we get $\frac{\epsilon \langle \phi \left(k\right)\rangle }{\langle k\rangle }\le \Theta \left(t\right)\le \frac{\langle \phi \left(k\right)\rangle }{\langle k\rangle }$. Then, from (4), we get

$$\frac{d{I}_k\left(t\right)}{dt}\le \lambda \left(k\right)\frac{\langle \phi \left(k\right)\rangle }{\langle k\rangle }(1-{\delta }_k)(N_k^{*}-I_k\left(t\right))-({\beta }_k+\gamma +d)I_k\left(t\right),t>0.$$

For any given constant $0<{\epsilon }_1<min\{\epsilon ,\frac{\lambda \left(k\right)(1-{\delta }_k)\frac{\langle \phi \left(k\right)\rangle }{\langle k\rangle }(1-N_k^{*})+{\beta }_k+\gamma +d}{\lambda \left(k\right)(1-{\delta }_k)\frac{\langle \phi \left(k\right)\rangle }{\langle k\rangle }+{\beta }_k+\gamma +d}\}$, there exists a $t_1>0$, when $t>t_1$, such that

$$I_k\left(t\right)\le A_k^1=\frac{\lambda \left(k\right)(1-{\delta }_k)\frac{\langle \phi \left(k\right)\rangle }{\langle k\rangle }{N}_k^{*}}{\lambda \left(k\right)(1-{\delta }_k)\frac{\langle \phi \left(k\right)\rangle }{\langle k\rangle }+{\beta }_k+\gamma +d}+{\epsilon }_1,t>t_1.$$

Then, from the second equation of (4), it implies that

$$\frac{d{Q}_k\left(t\right)}{dt}={\beta }_k(N_k^{*}-S_k\left(t\right)-Q_k\left(t\right))-(\eta +d)Q_k\left(t\right)\le {\beta }_k(N_k^{*}-Q_k\left(t\right))-(\eta +d)Q_k\left(t\right),t>0.$$

Hence, for any given constant $0<{\epsilon }_2<min\{\frac{1}{2},{\epsilon }_1,\frac{\eta +d+{\beta }_k-{\beta }_k{N}_k^{*}}{\eta +{\beta }_k}\}$, there exists a $t_2>t_1$, such that

$$Q_k\left(t\right)\le B_k^1=\frac{{\beta }_k{N}_k^{*}}{\eta +d+{\beta }_k},t>t_2.$$

On the other hand, from the first equation of (4), it can be inferred that

$$\frac{d{I}_k\left(t\right)}{dt}\ge \lambda \left(k\right)\frac{\epsilon \langle \phi \left(k\right)\rangle }{\langle k\rangle }(1-{\delta }_k)(N_k^{*}-I_k\left(t\right)-B_k^1)-({\beta }_k+\gamma +d)I_k\left(t\right).$$

Hence, for any given constant $0<{\epsilon }_3<min\{\frac{1}{3},{\epsilon }_2,\frac{\lambda \left(k\right)(1-{\delta }_k)\frac{\epsilon \langle \phi \left(k\right)\rangle }{\langle k\rangle }(N_k^{*}-B_k^1)}{\lambda \left(k\right)(1-{\delta }_k)\frac{\epsilon \langle \phi \left(k\right)\rangle }{\langle k\rangle }+{\beta }_k+\gamma +d}\}$, there exists a $t_3>t_2$, such that,

$$I_k\left(t\right)\ge a_k^1=\frac{\lambda \left(k\right)(1-{\delta }_k)\frac{\epsilon \langle \phi \left(k\right)\rangle }{\langle k\rangle }(N_k^{*}-B_k^1)}{\lambda \left(k\right)(1-{\delta }_k)\frac{\epsilon \langle \phi \left(k\right)\rangle }{\langle k\rangle }+{\beta }_k+\gamma +d}-{\epsilon }_3,t>t_3.$$

Next, following from the second equation of system (4), we also get

$$\frac{d{Q}_k\left(t\right)}{dt}={\beta }_k{I}_k\left(t\right)-(\eta +d)Q_k\left(t\right)\ge {\beta }_k{a}_k^1-(\eta +d)Q_k\left(t\right).$$

Similarly, for any given constant $0<{\epsilon }_4<min\{\frac{1}{4},{\epsilon }_3,\frac{{\beta }_k{a}_k^1}{\eta +d}\}$, there exists a $t_4>t_3,$ such that

$$Q_k\left(t\right)\ge b_k^1=\frac{{\beta }_k{a}_k^1}{\eta +d}-{\epsilon }_4,t>t_4.$$

Since ${\epsilon }_i$ is a small constant, it follows that

$$0<a_k^1\le I_k\left(t\right)\le A_k^1<1,0<b_k^1\le Q_k\left(t\right)\le B_k^1<1,t>t_4.$$

and, we also get $0<\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_k^1\le \Theta \left(t\right)\le \frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)A_k^1<\frac{\langle \phi \left(k\right)\rangle }{\langle k\rangle },t>t_4$. Again, by system (4), one has

$$\frac{d{I}_k\left(t\right)}{dt}\le \lambda \left(k\right)(1-{\delta }_k)(N_k^{*}-I_k\left(t\right)-b_k^1)\frac{1}{\langle k\rangle }\sum _{i=1}^n\phi \left(i\right)p\left(i\right)A_k^1-({\beta }_k+\gamma +d)I_k\left(t\right),t>t_4.$$

Hence, $0<{\epsilon }_5<min\{\frac{1}{5},{\epsilon }_4,\frac{(1-N_k^{*}+b_k^1)\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)A_k^1+{\beta }_k+\gamma +d}{\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)A_k^1+{\beta }_k+\gamma +d}\}$, there exists a $t_5>t_4,$ such that,

$$I_k\left(t\right)\le A_k^2=\frac{(N_k^{*}-b_k^1)\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)A_k^1}{\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)A_k^1+{\beta }_k+\gamma +d}+{\epsilon }_5,t>t_5.$$

Thus,

$$\frac{d{Q}_k\left(t\right)}{dt}={\beta }_k{I}_k\left(t\right)-(\eta +d)Q_k\left(t\right)\le {\beta }_k{A}_k^2-(\eta +d)Q_k\left(t\right),t>t_5$$

for any given constant $0<{\epsilon }_6<min\{\frac{1}{6},{\epsilon }_5,\frac{\eta +d-{\beta }_k{A}_k^2}{\eta +d}\}$, there exists a $t_6>t_5,$ such that,

$$Q_k\left(t\right)\le B_k^1=min\{B_k^1,\frac{{\beta }_k{A}_k^2}{\eta +d}+{\epsilon }_6\},t>t_6.$$

Consequently, one obtains that

$$\frac{d{I}_k\left(t\right)}{dt}\ge \lambda \left(k\right)(1-{\delta }_k)(N_k^{*}-I_k\left(t\right)-B_k^2)\frac{1}{\langle k\rangle }\sum _{i=1}^n\phi \left(i\right)p\left(i\right)a_k^1-({\beta }_k+\gamma +d)I_k\left(t\right),t>t_6.$$

Hence, for any given constant $0<{\epsilon }_7<min\{\frac{1}{7},{\epsilon }_6,\frac{(N_k^{*}-B_k^2)\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_k^1}{\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_k^1+{\beta }_k+\gamma +d}\}$, there exists a $t_7>t_6,$ such that,

$$I_k\left(t\right)\ge a_k^2=\frac{(N_k^{*}-B_k^2)\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_k^1}{\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_k^1+{\beta }_k+\gamma +d}-{\epsilon }_7,t>t_7.$$

Again, one has

$$\frac{d{Q}_k\left(t\right)}{dt}={\beta }_k{I}_k\left(t\right)-(\eta +d)Q_k\left(t\right)\ge {\beta }_k{a}_k^2-(\eta +d)Q_k\left(t\right),t>t_7.$$

for any given constant $0<{\epsilon }_8<min\{\frac{1}{8},{\epsilon }_7,\frac{{\beta }_k{a}_k^2}{\eta +d}\}$, there exists a $t_8>t_7,$ such that

$$Q_k\left(t\right)\ge b_k^2=max\{b_k^1,\frac{{\beta }_k{a}_k^2}{\eta +d}-{\epsilon }_8\},t>t_8.$$

It is clear that

$$0<a_k^1<a_k^2\le I_k\left(t\right)\le A_k^2<A_k^1<1,0<b_k^1<b_k^2\le Q_k\left(t\right)\le B_k^2<B_k^1<1,t>t_8.$$

Similarly, the calculation’s step m can be executed, resulting in four sequences. $\left\{A_k^m\right\},\left\{B_k^m\right\},\left\{a_k^m\right\},\left\{b_k^m\right\}$ are obtained:

$$\begin{array}{ccc}\hfill I_k\left(t\right)& \le & A_k^m=\frac{(N_k^{*}-b_k^{m-1})\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)A_k^{m-1}}{\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)A_i^{m-1}+({\beta }_k+\gamma +d)}+{\epsilon }_{4m-3},t>t_{4m-3},\hfill \\ \hfill Q_k\left(t\right)& \le & B_k^m=min\{B_k^{m-1},\frac{{\beta }_k{A}_i^m}{\eta +d}+{\epsilon }_{4m-2}\},t>t_{4m-2},\hfill \\ \hfill I_k\left(t\right)& \ge & a_k^m=\frac{(N_k^{*}-B_k^m)\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_i^{m-1}}{\lambda \left(k\right)(1-{\delta }_k)\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_i^{m-1}+({\beta }_k+\gamma +d)}-{\epsilon }_{4m-1},t>t_{4m-1},\hfill \\ \hfill Q_k\left(t\right)& \ge & b_k^m=max\{b_k^{m-1},\frac{{\beta }_k{a}_k^m}{\eta +d}-{\epsilon }_{4m}\},t>t_{4m}.\hfill \end{array}$$

It is clear that

$$0<a_k^m\le I_k\left(t\right)\le A_k^m<1,0<b_k^m\le Q_k\left(t\right)\le B_k^m<1,t>t_{4m}.$$

(13)

We get $\left\{a_k^m\right\},\left\{b_k^m\right\}$ are monotone increasing sequences and $\left\{A_k^m\right\},\left\{B_k^m\right\}$ are strictly monotone decreasing sequences. It is clear that

$$B_k^m=\frac{{\beta }_k{A}_k^m}{\eta +d}+{\epsilon }_{4m-2},b_k^m=\frac{{\beta }_k{a}_k^m}{\eta +d}-{\epsilon }_{4m}.$$

(14)

Since all of them are bounded monotonic sequences, the sequential limits exist. Let

$$\underset{t\to +\infty }{lim}{A}_k^{\left(m\right)}\left(t\right)=A_k,\underset{t\to +\infty }{lim}{B}_k^{\left(m\right)}\left(t\right)=B_k,\underset{t\to +\infty }{lim}{a}_k^{\left(m\right)}\left(t\right)=a_k,\underset{t\to +\infty }{lim}{b}_k^{\left(m\right)}\left(t\right)=b_k.$$

(15)

Since $\left\{{\epsilon }_m\right\}$ satisfies that $0<{\epsilon }_m<\frac{1}{m}$, one has ${\epsilon }_m\to 0$ as $m\to \infty$. Therefore, from (14), as $m\to \infty$, we get

$$A_k=\frac{(N_k^{*}-b_k)\lambda \left(k\right)(1-{\delta }_k)A}{\lambda \left(k\right)(1-{\delta }_k)A+({\beta }_k+\gamma +d)},B_k=\frac{{\beta }_k{A}_k}{\eta +d},$$

$$a_k=\frac{(N_k^{*}-B_k)\lambda \left(k\right)(1-{\delta }_k)a}{\lambda \left(k\right)(1-{\delta }_k)a+({\beta }_k+\gamma +d)},b_k=\frac{{\beta }_k{a}_k}{\eta +d}.$$

where $A=\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)A_i,a=\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_i$ Then, we get

$$\begin{array}{cc}\hfill A_k=& \frac{N_k^{*}{\varsigma }_{k}A(\eta +d)[a{\varsigma }_k(\eta +d-{\beta }_k)+(\eta +d)({\beta }_k+\gamma +d)]}{H_k},\hfill \\ \hfill a_k=& \frac{N_k^{*}{\varsigma }_{k}a(\eta +d)[A{\varsigma }_k(\eta +d-{\beta }_k)+(\eta +d)({\beta }_k+\gamma +d)]}{H_k}.\hfill \end{array}$$

(16)

where ${\varsigma }_k=\lambda \left(k\right)(1-{\delta }_k)$,$H_k={(\eta +d)}^2({\varsigma }_{k}A+({\beta }_k+\gamma +d))({\varsigma }_{k}a+({\beta }_k+\gamma +d))-{\varsigma }_k^2{\beta }_k^{2}Aa.$ For any k, if $\eta +d={\beta }_k$, it is clear that $A=a$. If $\eta +d\ne {\beta }_k$, Substituting (16) into A:

$$\begin{array}{cc}\hfill A=& \frac{1}{\langle k\rangle }\sum _{i=1}^n\phi \left(i\right)p\left(i\right)\frac{N_i^{*}{\varsigma }_{i}A(\eta +d)[a{\varsigma }_i(\eta +d-{\beta }_i)+(\eta +d)({\beta }_i+\gamma +d)]}{H_i},\hfill \\ \hfill a=& \frac{1}{\langle k\rangle }\sum _{i=1}^n\phi \left(i\right)p\left(i\right)\frac{N_i^{*}{\varsigma }_{i}a(\eta +d)[A{\varsigma }_i(\eta +d-{\beta }_i)+(\eta +d)({\beta }_i+\gamma +d)]}{H_i}.\hfill \end{array}$$

(17)

We get $\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)\frac{{\varsigma }_i{N}_i^{*}}{H_i}\left[{(\eta +d)}^2({\beta }_i+\gamma +d)(A-a)\right]=0$. Therefore, $A=a$. Above all, from $A=a$, which is equivalent to $A_i=a_i$ for any $i=1,2,\dots ,n$. Then, from (13) and (14), it follows that

$$I_k^{*}=\underset{t\to +\infty }{lim}{I}_k\left(t\right)=A_k=a_k=\frac{(N_k^{*}-b_k)\lambda \left(k\right)(1-{\delta }_k){\Theta }^{*}}{\lambda \left(k\right)(1-{\delta }_k){\Theta }^{*}+({\beta }_k+\gamma +d)},$$

$$Q_k^{*}=\underset{t\to +\infty }{lim}{Q}_k\left(t\right)=B_k=b_k=\frac{{\beta }_k{I}_k^{*}}{\eta +d}.$$

which is consistent with Theorem 2. This completes the proof. □

5. The Optimal Control for SIQS Model

In this section, we study the optimal control for system (1). To meet the control objective and minimize control expenses at the same time, a viable approach to target epidemic outbreaks is by implementing time-varying control using optimal control theory for various infectious diseases, for example, the results in Culshaw et al. [44], Chen et al. [45] and Abboubakar et al. [46]. The aim of this section is to identify the optimal quarantine measure for controlling the transmission of diseases.

First, we introduce an time-varying control rate $u_k\left(t\right)$ represents the percentage of infected individuals with degree k being quarantined. Hence, system (3) can be written as

$$\left\{\begin{array}{c}\frac{d{S}_k\left(t\right)}{dt}=bk[1-N_k^{*}]{\Psi }^{*}-d{S}_k\left(t\right)-\lambda \left(k\right)(1-{\delta }_k)S_k\left(t\right)\Theta \left(t\right)+\gamma I_k\left(t\right)+\eta Q_k\left(t\right),\hfill \\ \frac{d{I}_k\left(t\right)}{dt}=\lambda \left(k\right)(1-{\delta }_k)S_k\left(t\right)\Theta \left(t\right)-(u_k\left(t\right)+\gamma +d)I_k\left(t\right),\hfill \\ \frac{d{Q}_k\left(t\right)}{dt}=u_k\left(t\right)I_k\left(t\right)-(\eta +d)Q_k\left(t\right).\hfill \end{array}\right.$$

(18)

We note the quarantine function $u\left(t\right)=(u_1\left(t\right),u_2\left(t\right),...u_n\left(t\right))$ are bounded, Lebesgue integrable functions. Since our aim is to reduce both the scale of infection and control strengths, we consider the following objective functional $J\left(u_k\right)$ associated with model in the control set $\Omega =\{(u_1\left(t\right),u_2\left(t\right),...u_n\left(t\right))\mid 0\le u_i\le 0.9\}$, $B_k$ denotes a positive weight parameter. The objective (cost) functional given by

$$J\left(u_k\right)={\int }_0^T\sum _{k=1}^n[I_k\left(t\right)+\frac{1}{2}{B}_k{u}_k^2\left(t\right)]\mathrm{d}t.$$

(19)

The Lagrange function for optimal control problems (18) and (19) is defined as

$$L=\sum _{i=1}^n[I_i+\frac{B_i}{2}{u}_i^2\left(t\right)].$$

In order to solve the optimal control problem, the existence of an optimal control must be assured. We firstly establish the existence of a solution for the control system (18) by applying the method in Refs. [45,46].

Theorem 5.

For the objective functional $J\left(u_k\right)$ associated with model (18) defined in Ω, there exists an optimal quarantine strategy ${\mu }^{*}=({\mu }_1,\dots ,{\mu }_n)$ minimizing $J\left(u\right)$.

Proof. 

(i) Clearly, for any $0\le t\le T$, the $\Omega$ is an nonempty set of Lebesgue integrable functions. (ii) It is obviously that the solutions are bounded ensuring the boundedness and convexity of admissible control set. (iii) Model (18) can be written as

$$\frac{d\Phi \left(t\right)}{dt}=F(t,\Phi ,u_k\left(t\right))$$

(20)

where $F(t,\Phi ,u_k)=(F_1,F_2,\dots ,F_n)$ is the right side of system (18), and $\Phi ={({\Phi }_1,{\Phi }_2,\dots ,{\Phi }_n)}^T$, ${\Phi }_k={(S_k\left(t\right),I_k\left(t\right),Q_k\left(t\right))}^T$. It follows that

$$\begin{array}{ccc}\hfill |F(t,{\Phi }_1,u_k\left(t\right))-F(t,{\Phi }_2,u_k\left(t\right))|& \le & \sum _i|F_i(t,{\Phi }_1,u_k\left(t\right))-F_i(t,{\Phi }_2,u_k\left(t\right))|.\hfill \end{array}$$

For any $k\in {\mathbb{N}}_n$, since $|{\Theta }_1\left(t\right)-{\Theta }_2\left(t\right)|=\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)|I_{i1}\left(t\right)-I_{i2}\left(t\right)|$ and the boundedness of solution, there exists a constant $M_1$ satisfied that $|{\Theta }_1\left(t\right)-{\Theta }_2\left(t\right)|\le M_1|I_{k1}\left(t\right)-I_{k2}\left(t\right)|$. Then, There exists a positive constant M to guarantee the inequality $|F_k\left({\Phi }_1\right)-F_k\left({\Phi }_2\right)|\le M|{\Phi }_{k1}-{\Phi }_{k2}|$ is satisfied. (iv) Let $0<\xi <1$ and $m_k\left(t\right),v_k\left(t\right)\in \Omega$, we get $L(t,\Phi ,(1-\xi )m\left(t\right)+\xi v\left(t\right))-(1-\xi )L(t,\Phi ,m\left(t\right))-\xi L(t,\Phi ,v\left(t\right))={\sum }_{i=1}^n\frac{B_i}{2}\xi (\xi -1){(v_i-m_i)}^2<0$. (v) There exist $c1>0$, $\rho >1$, $c2>0$ such that $L(t,\Phi ,u_k\left(t\right))={\sum }_{i=1}^n[I_i+\frac{B_i}{2}{u}_i^2\left(t\right)]\ge {\sum }_{i=1}^n\frac{B_i}{2}{u}_i^2\left(t\right)\ge \frac{B_{min}}{2}{\sum }_{i=1}^n{u}_i^2\left(t\right)\ge C_1{\left|u\right|}^2-C_2$.

By utilizing the results of Refs. [45,46], we get that there exists a optimal control ${\mu }^{*}=({\mu }_1,\cdots ,{\mu }_n)$ minimizing $J\left(\mu \right)$. The proof is completed. □

The solution to the optimal control problem will be determined through Pontryagin’s minimum principle [47]. We note the Hamiltonian H as

$$H=\sum _{i=1}^n[I_k\left(t\right)+\frac{1}{2}{B}_k{u}_k^2\left(t\right)]\mathrm{d}t+\sum _{i=1}^3\sum _{k=1}^n{\omega }_{ki}{F}_{ki}(t,\Phi ,u_k\left(t\right)).$$

(21)

If $({\Phi }^{*},u_k^{*})$ is an optimal solution of the optimal control problem, then there exists a non-trivial vector function ${\omega }_k=({\omega }_{k1},{\omega }_{k2},{\omega }_{k3})$, $k=1,2,\dots ,n$ satisfying the following equalities:

$$\left\{\begin{array}{c}\frac{d\Phi \left(t\right)}{dt}=\frac{\partial H(t,{\Phi }^{*}\left(t\right),u_k^{*},\omega \left(t\right))}{\partial \omega },\hfill \\ \frac{\partial H(t,{\Phi }^{*}\left(t\right),u_k^{*},\omega \left(t\right))}{\partial u_k\left(t\right)}=0,\hfill \\ \frac{d\omega \left(t\right)}{dt}=-\frac{\partial H(t,{\Phi }^{*}\left(t\right),u_k^{*},\omega \left(t\right))}{\partial \Phi }.\hfill \end{array}\right.$$

(22)

It follows from the derivation above

$$\left\{\begin{array}{c}{\omega }^{*}=0,\frac{\partial H}{\partial u_k\left(t\right)}<0,\hfill \\ 0<{\omega }^{*}<0.9,\frac{\partial H}{\partial u_k\left(t\right)}=0,\hfill \\ {\omega }^{*}=0.9,\frac{\partial H}{\partial u_k\left(t\right)}>0.\hfill \end{array}\right.$$

(23)

Now the necessary conditions are implemented on the Hamiltonian H in (21) and we obtain the following results.

Theorem 6.

For a given optimal control solution $\mu \left(t\right)$ and the corresponding system solution $S_k^{*}\left(t\right),Q_k^{*}\left(t\right)$, $I_k^{*}\left(t\right)$, there is an adjacency function ${\omega }_{k1}\left(t\right),{\omega }_{k2}\left(t\right),{\omega }_{k3}\left(t\right)$ that satisfies the following conditions

$$\left\{\begin{array}{cc}\hfill \frac{d{\omega }_{k1}\left(t\right)}{dt}& =\lambda \left(k\right)(1-{\delta }_k){\Theta }^{*}\left(t\right)[{\omega }_{k1}\left(t\right)-{\omega }_{k2}\left(t\right)]+d{\omega }_{k1}\left(t\right)-1,\hfill \\ \hfill \frac{d{\omega }_{k2}\left(t\right)}{dt}& =\frac{\lambda \left(k\right)(1-{\delta }_k)S_k^{*}\left(t\right)\phi \left(k\right)p\left(k\right)}{\langle k\rangle }({\omega }_{k1}\left(t\right)-{\omega }_{k2}\left(t\right))-\gamma {\omega }_{k1}\left(t\right)\hfill \\ \hfill & +(u_k+\gamma +d){\omega }_{k2}\left(t\right)-u_k{\omega }_{k3}\left(t\right)-1,\hfill \\ \hfill \frac{d{\omega }_{k3}\left(t\right)}{dt}& =-\eta {\omega }_{k1}\left(t\right)+(\eta +d){\omega }_{k3}\left(t\right).\hfill \end{array}\right.$$

(24)

with conditions ${\omega }_{k1}\left(T\right)$=${\omega }_{k2}\left(T\right)$=${\omega }_{k3}\left(T\right)$=0. Furthermore, the optimal control $u_k^{*}=(u_1^{*},u_2^{*},...,u_n^{*})$ is given by

$$\begin{array}{cc}\hfill u_k^{*}=& min\{max[0,\frac{I_k^{*}\left(t\right)({\omega }_{k2}-{\omega }_{k3})}{B_k}],0.9\}\hfill \end{array}$$

(25)

Proof. 

We apply for the necessary conditions in the Pontryagin maximum principle with Hamiltonian function (21) to derive the adjoint variables is determined by the following equations.

$$\begin{array}{ccc}\hfill \frac{d{\omega }_{k1}\left(t\right)}{dt}& =& \lambda \left(k\right)(1-{\delta }_k){\Theta }^{*}\left(t\right)[{\omega }_{k1}\left(t\right)-{\omega }_{k2}\left(t\right)]+d{\omega }_{k1}\left(t\right)-1,\hfill \\ \hfill \frac{d{\omega }_{k2}\left(t\right)}{dt}& =& \frac{\lambda \left(k\right)(1-{\delta }_k)S_k^{*}\left(t\right)\phi \left(k\right)p\left(k\right)}{\langle k\rangle }({\omega }_{k1}\left(t\right)-{\omega }_{k2}\left(t\right))-\gamma {\omega }_{k1}\left(t\right)\hfill \\ & & +(u_k+\gamma +d){\omega }_{k2}\left(t\right)-u_k{\omega }_{k3}\left(t\right)-1,\hfill \\ \hfill \frac{d{\omega }_{k3}\left(t\right)}{dt}& =& -\eta {\omega }_{k1}\left(t\right)+(\eta +d){\omega }_{k3}\left(t\right).\hfill \end{array}$$

By the optimal conditions, we have

$$\frac{\partial H}{\partial u_k\left(t\right)}=B_k{u}_k^{*}-{\omega }_{k2}{I}_k^{*}+{\omega }_{k3}{I}_k^{*}=0.$$

It follows that ${\mu }_k^{*}=\frac{({\omega }_{k2}-{\omega }_{k3})I_k^{*}}{B_k}$. Using the property of the control space, we obtain

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\omega }^{*}=0,\frac{({\omega }_{k2}-{\omega }_{k3})I_k^{*}}{B_k}<0,\hfill \\ 0<{\omega }^{*}<0.9,\frac{({\omega }_{k2}-{\omega }_{k3})I_k^{*}}{B_k}=0,\hfill \\ {\omega }^{*}=0.9,\frac{({\omega }_{k2}-{\omega }_{k3})I_k^{*}}{B_k}>0.\hfill \end{array}\right.\end{array}$$

Then, we have the optimal control $u_k^{*}$ in compact notation,

$$\begin{array}{c}\hfill u_k^{*}=min\{max[0,\frac{({\omega }_{k2}-{\omega }_{k3})I_k^{*}}{B_k}],0.9\}\end{array}$$

That completes the proof. □

6. Simulations

In this section, we conducted some numerical simulations to validate and supplement our theoretical results, and to study the impact of parameters on the dynamics of disease, in order to find better control strategies. The numerical simulations include three parts, the first part is to verify the main theoretical results including the epidemic threshold and global stability. We will discuss and make some simulations to show the effectiveness of different immunity and quarantine strategies in second part. The numerical simulation in the third part demonstrates the effectiveness of the optimal control method that has been proved in Section 5. We constructed a BA (Barab$\acute{a}$ si-Albert) scale-free network [48], which satisfies a power-law degree distribution $P\left(k\right)\sim C{k}^{-3}$, the network size N is set as 1000, $n=k_{max}$=80. This network evolved from the initial network with size of $m_0=3$ and during each time step, new nodes with $m=3$ edges are added to the network. Without loss of generality, we note $\lambda \left(k\right)=\lambda k$.

(I) Firstly, we perform some simulations in Figure 2 to verify the important results of global stability of equilibria. We denote

$$\begin{array}{c}\hfill \phi \left(k\right)=\left\{\begin{array}{c}k,k\le 30;\hfill \\ 30,k>30.\hfill \end{array}\right.\end{array}$$

Figure 2(a) shows the relationship between $I\left(t\right)$ and infectious rate $\lambda$ on some different networks. Based on the basic reproduction number $R_0$, we can acquire the epidemic threshold of the disease is ${\lambda }_c=\langle k\rangle \setminus {\sum }_{i=1}^n\frac{\phi \left(i\right)p\left(i\right)i{N}_i^{*}(1-{\delta }_i)}{{\beta }_i+\gamma +d}$, where $b=0.2$, $d=0.16$, ${\beta }_k=0.17$, $\eta =0.3$, ${\delta }_k=0.2$, $\gamma =0.36$. One can see that the epidemic threshold is consistent with our theoretical results very well. Figure 2(b) depicts that the disease-free equilibrium cannot undergo a Hopf bifurcation, i.e., if $R_0<1$, the disease free equilibrium is global stability, and if $R_0>1$, then the endemic equilibrium is persist and global stability, where $b=0.2$, $d=0.1$, ${\beta }_k=0.2$, $\eta =0.3$, ${\delta }_k=0.2$, $\lambda$=0.01, 0.038, 0.055, 0.21, 0.28, 0.35 and $\gamma$=0.36, 0.36, 0.36, 0.18, 0.18, 0.18.

(II) Secondly, we depicts the effectiveness of the immune and quarantining control strategies in Figure 3 to Figure 4. In this part, we also discuss the SIQS model (1) on the network with various immunization (proportional immune strategy $\tilde{\delta }$, target immune strategy ${\delta }_k^t$, acquaintance immune strategy ${\delta }_k^{ac}$) schemes based on [49] and define quarantine schemes of infected individual (proportional quarantine strategy ${\beta }_k$, target quarantine strategy ${\beta }_k^t$, acquaintance quarantining strategy ${\beta }_k^{ac}$) according to various immunization schemes. Then, we define the heterogeneous quarantine rate ${\beta }_k$ as follows,

Proportional quarantine strategy: Denote the average quarantine rate of proportional quarantine $\tilde{\beta }$, $0<\tilde{\beta }<1$. Randomly selecting one infected node for isolation, the quarantine rate is independent of the degree of node, which is also a situation discussed in many papers.

Target quarantine strategy: we can devise a quarantine strategy for the infected nodes according to the definition of immunization [49]. Introduce an upper threshold $\kappa$, such that all infected nodes with connectivity $k>\kappa$ are prioritize quarantined, i.e., we defined the quarantine rate ${\beta }_k^t$ by

$$\begin{array}{c}\hfill {\beta }_k^t=\left\{\begin{array}{c}1,k>\kappa ;\hfill \\ c,k=\kappa ;\hfill \\ 0,k<\kappa .\hfill \end{array}\right.\end{array}$$

where $0<c\le 1$, and ${\sum }_k{\beta }_k^{t}P\left(k\right)=\widetilde{{\beta }_k^t}$, where $\widetilde{{\beta }_k^t}$ is the average quarantine rate of target quarantine.

Acquaintance quarantining strategy: Select a random portion p from the N nodes. The likelihood of quarantining a infected node with degree k is given by $\frac{kP\left(k\right)}{N\langle k\rangle }$, therefore ${\beta }_k^{ac}=pkP\left(k\right)/\langle k\rangle$. We note that $\widetilde{{\beta }_k^{ac}}$ denote the average quarantine rate of acquaintance quarantine, where $\widetilde{{\beta }_k^{ac}}={\sum }_k{\beta }_k^{ac}P\left(k\right)$.

Figure 3(a) shows the $I\left(t\right)$ under different quarantining control strategies, the target quarantine strategy is the most effective in controlling disease. Figure 3(b) compares the effectiveness of different quarantine strategies. We get that the target quarantine strategy has better effectiveness than others. Figure 3(c) shows the effectiveness of different immunity strategies, it shows that all three immunization schemes are effective compared to the case without any immunization; the targeted immunization scheme is more efficient than the proportional scheme discussed.

Figure 4(a) depicts the average infectious density with respect to quarantine rate under different immunity strategies. It also shows that the target strategy is better than proportional and acquaintance. We further show the average infectious density with respect to infected rate for the immunization schemes and quarantine strategies targeted, and acquaintance in Figure 4(b). We get that the targeted immunity and the target quarantine strategy have better effectiveness than others.

(III) Thirdly, we analyze and illustrate the optimal control and various quarantine schemes by numerical simulations in Figure 5 to Figure 6. In Figure 5, we depict $I\left(t\right)$ under different control strategies. In order to provide a clearer representation of the findings in Figure 5(a), we additionally computed the objective function $J\left(\mu \right)$ values for diverse control approaches, as shown in Table 2. In the absence of any control strategy (i.e., $u_k=0$), the infection ultimately breaks out and reaches a stable level of infection. We also set a fixed control strategy compared with the optimal control strategy, which can find that if we want to obtain better results than the optimal control strategy, it will usually lead to a significant increase in control costs (i.e., $u_k^{max}=0.9$). In Figure 5(b), we demonstrate the comparative effects of optimal control, without control, and acquaintance immunity combining with multiple quarantine strategies on infectious disease control. Compared with other strategies, the optimal control strategy has the best control effect.

Figure 6(a) illustrates the $I\left(t\right)$ under optimal control under three cases: Low-cost at $B_k=0.01$, moderate-cost at $B_k=1$, and high-cost at $B_k=100$. Even under high costs, as seen in the third scenario, the scale of infection can be notably curtailed through optimal control. Figure 6(b) displays the fluctuation of ${\mu }_k$ across varying degrees under optimal regulation. It is evident that ${\mu }_k$ tends to escalate as the degree rises. This indicates a higher level of infection, resulting in a higher value of ${\mu }_k$.

7. Conclusion

We propose an SIQS epidemic model with heterogeneous immune and quarantine measures on heterogeneous networks. We have performed a detailed mathematical analysis of our system, which deduces that the basic reproduction number $R_0$, is collectively determined by the infection rate and the effectiveness of disease control measures. Specifically, when $R_0<1$, the disease-free equilibrium $E^0$ exhibits global asymptotic stability. When $R_0>1$, the endemic equilibrium $E^{*}$ demonstrates both local and global asymptotic stability. At the same time, we conducted an in-depth theoretical examination on system (1) to tackle the optimal quarantine control issue, confirming the presence of optimal solutions. We define different forms of quarantine measures and numerically simulated the effectiveness of a combination of different quarantine strategies and immune schemes in controlling diseases. Our findings indicate that the implementation of any control strategy can have a favorable impact on the magnitude of an epidemic outbreak, with the optimal control being particularly effective. Among them, the effect of target quarantine scheme can be comparable to the effect of optimal control. However, target isolation requires knowing the degree of all nodes in the network, which is a significant cost and workload.

A variety of control strategies, such as quarantine and immune schemes, can be employed for the effective management, mitigation, and potential eradication of infectious diseases. We expect that this paper can contributes to a understanding of the dynamics of infectious disease with the control measures. However, our paper is not based on a specific infectious disease, so there are still certain limitations in providing a specific basis for public health prevention measures. Studying the control measures for a specific disease based on infectious disease data will be a very meaningful work.