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

Abstract

Complex network theory involves network structure and dynamics; dynamics on networks and interactions between networks; and dynamics developed over a network. As a typical application of complex networks, the dynamics of disease spreading and control strategies on networks have attracted widespread attention from researchers. We investigate the dynamics and optimal control for an epidemic model with demographics and heterogeneous asymmetric control strategies (immunization and quarantine) on complex networks. We derive the epidemic threshold and study the global stability of disease-free and endemic equilibria based on different methods. The results show that the disease-free equilibrium cannot undergo a Hopf bifurcation. We further study the optimal control strategy for the complex system and obtain its existence and uniqueness. Numerical simulations are conducted on scale-free networks to validate and supplement the theoretical results. The numerical results indicate that the asymmetric control strategies regarding time and degree of node for populations are superior to symmetric control strategies when considering control cost and the effectiveness of controlling infectious diseases. Meanwhile, the advantages of the optimal control strategy through comparisons with various baseline immunization and quarantine schemes are also shown.

1. Introduction

The prevalence and outbreak of infectious diseases are always important issues affecting the national economy and people’s livelihoods. Research on the transmission mechanism and dynamics of infectious diseases has 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 [1,2]. 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 [3,4,5,6]. Based on different mechanisms of disease transmission, Kermack and McKendrick introduced the SIS and SIR epidemic models [7,8,9]. 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 epidemic transmission dynamics and control measures may have a significant impact on the prevention and control of infectious diseases, and this topic has attracted widespread attention in scientific research [10,11,12]. During the spread of infectious diseases, the infected transmission mainly occurs through contact between susceptible and infected individuals. Therefore, quarantine, which prevents this contact, has become a commonly employed strategy for controlling disease. In the past two decades, some important methods [13,14,15,16] have been introduced to control outbreaks of infectious diseases on the assumption of a uniformly mixing population, which show that the quarantine rate might be a factor contributing to the observed sustained oscillations in some directly transmitted viral diseases. Since Barabási and Albert [17] proposed a scale-free network model, in which the degree distribution follows a power-law distribution, researchers have gradually focused on the study of transmission and control strategies of infectious diseases on complex networks and obtained a lot of useful and insightful results. The most influential results are the SIS and SIR epidemic models established by Pastor Satoras and Vespigani [18] through mean field theory, which indicate that the epidemic threshold of a disease will tend to 0 when the network scale is large enough. Furthermore, there are many results in theoretical research on disease control strategies on heterogeneous networks. Many researchers have introduced quarantine compartments Q into some epidemic models (e.g., SIS model [10,19,20], SIR model [10], SIRS model [21], SEIR model [22,23], etc.) on scale-free networks; these results all suggest that the epidemic threshold of a disease is closely linked to the network’s structural characteristics and quarantine rate, and the heterogeneity of networks and higher infectivity increase the risk of disease transmission, while quarantine measures help to prevent infection. Therefore, constructing a network-based transmission model can more realistically reflect the propagation laws of infectious diseases in real contact networks and explore the coupling effects between network structure, epidemic dynamics, and control strategies of diseases. However, those models either ignore the demographics or assume that infectious disease control strategies may maintain symmetry among individuals, that is, the control strategy for the nodes in the network is uniform.

Control measures for infectious diseases may maintain symmetry in time by taking the same measure with the same strength at different time points. This can include isolation, social distancing, or other control measures. Such strategies can effectively reduce the scale of infection. However, relying solely on the best uniform control measures feasible given economic, social, or other constraints may not be sufficient to achieve the disease control goals. An attractive alternative approach is the use of optimal control theory to employ time-varying control within certain bounds to strike a balance between control objectives and the associated costs. Optimal control applications in epidemic dynamics have primarily concentrated on homogeneous contact networks [24,25,26,27,28,29]. Recently, many intriguing works on the optimal control of epidemics on heterogeneous networks have emerged [21,30,31,32,33,34,35,36,37]. Li et al. presented a nonlinear SIQS epidemic model of networks and explored the issue of optimal quarantine control to minimize the cost of control measures [34]. Zhang et al. examined the optimal control of an SIQRS model that includes vaccination in a network and studied the effects of different control strategies [36]. Yang et al. investigated the stability and optimal control of SIS epidemic systems in directed networks [37]. However, these contributions were restricted to an optimal control strategy for isolation and did not compare results 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. Based on the above considerations, we study an SIS epidemic model with heterogeneous immunization for susceptible individuals and quarantine for infected individuals, while also considering the demographics of individuals on heterogeneous networks. The organization of this paper is as follows. In Section 2, we introduce an SIQS model with heterogeneous asymmetric control strategies. Traditional infectious disease control strategies maintain the symmetry principle that the quarantine rate of an infected node is constant across individuals. By heterogeneous asymmetric control strategies, we mean that the quarantine and immunization rates of an infected node are not constant but are related to the node degree. In other words, the control strategy on the network is not uniform, and instead, different control strategies are formulated based on degree, on complex networks, and in this section, we state our assumptions. In Section 3, we obtain the equilibria and basic reproduction number. The main results are shown in Section 4, where the uniform persistence and global stability of the disease-free and endemic equilibria of a system are analyzed by different mathematical methods. The optimal time-varying quarantine control for the model is considered in Section 5. In Section 6, some numerical simulations are conducted to illustrate and supplement the analysis results. Finally, a brief conclusion is given in Section 7.

2. Description and Formation of Epidemic Models

Guided by the symmetry of the contact mode on real contact networks when diseases spread, we consider a population situated on an undirected complex network N. Every node has three optional states: the susceptible state (S), quarantined state (Q), and infected state (I) [20]. We divided the population into n groups. Let $S_k\left(t\right)$, $I_k\left(t\right)$, and $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. A schematic diagram is shown in Figure 1, and the meanings of the parameters are listed in

Table 1. Symbols employed in models.

Symbols	Description
$p\left(k\right)$	Proportion of nodes with degree k.
$\langle k\rangle$	Average degree  $\left(\langle k\rangle ={\sum }_{k}kp\left(k\right)\right)$.
n	Maximum degree.
b	Birth rate.
d	Natural death rate.
${\Psi }_k\left(t\right)$	Fertile contact probability between a node with degree k and its neighbors.
$\lambda \left(k\right)$	Transmission rate of infected nodes with degree k.
${\delta }_k$	Vaccination rate of susceptible nodes with degree k.
${\beta }_k$	Quarantine rate of infected nodes with degree k.
$\gamma$	Recovery rate of infected nodes.
$\eta$	Recovery rate of quarantined nodes.

. These factors are as follows:

(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 nonvacant node, the vacant node i will be activated to generate a new susceptible node at the birth rate b. Each nonvacant node becomes a vacant node at a natural death rate d per unit time. We assume that each nonvacant 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 at the immune rate ${\delta }_k$. The infected nodes with degree k will be quarantined at rate ${\beta }_k$. The quarantined individuals will recover to a susceptible node at rate $\eta$. Nodes with the same degree have identical quarantine and immunization strategies, while those with different degrees have different strategies.

(3) 

SIS epidemic framework: Infection $S\to I$: At the initial moment, some nodes are randomly selected as infected nodes. At each time step, the possibility that each infected node i will 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$ [38,39,40], $\phi \left(k\right)$ = A [41], $\phi \left(k\right)$ = $k^m$ [42], $\phi \left(k\right)$ = $\frac{a{k}^m}{1+b{k}^m}$ [43]. 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 a 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 the above assumptions, 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, …, 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 a link which originates from a node with degree k points to an infected node. ${\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 possibility of selecting a neighbor with degree i of an empty node to activate it. On uncorrelated networks, the likelihood is 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)$, and $Q\left(t\right)={\sum }_{k=1}^n{Q}_k\left(t\right)p\left(k\right)$.

Let $N_k\left(t\right)$ denote 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\left(t\right)$ 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. [44], we obtain ${lim}_{t\to \infty }{N}_k\left(t\right)=0$ when $b\le d$, and 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 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$ is a solution to 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$, and $k\in {\mathbb{N}}_n$ for any $t\ge 0$, i.e., Γ is positively invariant.

Proof. 

∀$x_0\in \Gamma$. Thus, $\forall t>0$, and so 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 find 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 Equation (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 exists a $j\in \{1,2,\dots ,n\}$ and $t_0\ge 0$, satisfying

$$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$; therefore, 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 $\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 obtain $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 [19]. 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)

Let X denote the space of functions and $x=(I_1\left(t\right),Q_1\left(t\right),\dots ,I_n\left(t\right),Q_n\left(t\right))\in X$ denote 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)

Proposition 1 can show that $\Gamma$ is positively invariant for system (4). All viable steady states of system (4) satisfy 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 a disease-free equilibrium of (4) $E^0=(0,0,\dots ,0))$, and the positive equilibrium $E^{*}={(I_k^{*},Q_k^{*})}_k$, where $E^{*}$ satisfies

$$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 satisfies

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

Equation (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$, we obtain

$$\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 $f^{\prime }(\Theta )<0$ and $\underset{\Theta \to +\infty }{lim}f(\Theta )=0$, it is clear that Equation (9) has a 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$, Equation (6) has a unique positive solution ${\Theta }^{*}$, namely, a unique positive equilibrium point $E^{*}$ of system (4) exists as described by Equation (7) if and only if $R_0>1$. The proof is completed.  □

Remark 1. 

Without losing generality, we know that there always exists constant ${\beta }_{min}$, ${\beta }_{max}$, ${\delta }_{min}$, ${\delta }_{max}$, $N_{min}$ $\in (0,1)$ that satisfy

${\beta }_{min}\le {\beta }_k\le {\beta }_{max},{\delta }_{min}\le {\delta }_k\le {\delta }_{max},N_{min}\le N_k^{*}<1.$ Then, let $\phi \left(k\right)=k,$ $\lambda \left(k\right)=\lambda k$, we can obtain $R_0\ge \lambda \frac{\langle k^2\rangle }{\langle k\rangle }\frac{(1-{\delta }_{max})N_{min}}{{\beta }_{max}+\gamma +d}.$ Based on [18], we can find when the network size is large enough and the second-order moment of degree $\langle k^2\rangle \to \infty$; then, $R_0$ can always be guaranteed to be greater than 1 and the disease will persist. Therefore, the topology of the network, such as the average degree $\langle k\rangle$ and the second moment of degree $\langle k^2\rangle$, has an important impact on the epidemic threshold of infectious diseases.

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 used in [45] 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 of order n, $A_{21}={\beta }_k{E}_n$, $A_{22}=-(\eta +d)E_n$, $E_n$ is the 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).$$

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 are determined by $|\rho E_n-A_{11}|=0$. Through some elementary transformations, $|\rho E_n-A_{11}|=0$ can be changed 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 obtain

$$(\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 that $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, 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, 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$, $i=2,\dots ,n$.

Furthermore, because $N\left(y\right)\le 0$, we have $\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 obtain

$$\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}$$

From Proposition 1, 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$. 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 obtain $\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 disease, which is important in proving the global stability of endemic equilibrium $E^{*}$. We obtain the uniform persistence based on Theorem 4.2 of Ref. [46]. The result of uniform persistence is shown 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 Theorem 4.2 of [46], 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 a positive eigenvector $\tilde{h}=(h_1,h_2,\dots ,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$, 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 the 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. All conditions are satisfied, and the proof of uniform persistence is completed.  □

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

Theorem 4. 

Suppose that $(I_k\left(t\right),Q_k\left(t\right))$ is a solution of system (4), satisfying the 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 globally attractive.

Proof. 

Firstly, from Lemma 1 and for any $k\in \{1,2,\dots ,n\}$, we find that there exists $0<\epsilon \ll 1,$ such that, for t that is sufficiently large, $I_k\ge \epsilon$ and $Q_k\ge \epsilon$ are satisfied. Hence, when $t>0$, we obtain $\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 system (4), we obtain

$$\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 system (4), we have

$$\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 system (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 obtain

$$\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 obtain $0<\frac{1}{\langle k\rangle }{\sum }_{i=1}^n\phi \left(i\right)p\left(i\right)a_i^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_i^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_i^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_i^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_i^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_i^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^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$$

and 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_i^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_i^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^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_i^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^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 as follows:

$$\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_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-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 find that $\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 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 obtain

$$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 obtain

$$\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, we obtain

$$\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 obtain $\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, $A=a$ 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 the 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 targeting epidemic outbreaks is by implementing time-varying control using optimal control theory for various infectious diseases, for example, the results in Culshaw et al. [47], Chen et al. [48], and Abboubakar et al. [49]. The aim of this section is to identify the optimal quarantine measures for controlling the transmission of diseases.

First, we introduce a time-varying control rate $u_k\left(t\right)$ representing 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 functions $u\left(t\right)=(u_1\left(t\right),u_2\left(t\right),\dots u_n\left(t\right))$ are bounded, Lebesgue integrable functions. Since our aim is both to reduce the scale of infection and to control its strength, we consider the following objective functional $J\left(u_k\right)$ associated with the model in the control set $\Omega =\{(u_1\left(t\right),u_2\left(t\right),\dots u_n\left(t\right))\mid 0\le u_i\le 0.9\}$. $B_k$ denotes a positive weight parameter. The objective (cost) functional is 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 applied Corollary 4.1 in Ref. [50] to prove the existence of the solution, which is shown in the proof of Theorem 5.

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 a nonempty set of Lebesgue integrable functions. (ii) It is obvious that the solutions are bounded, ensuring the boundedness and convexity of the 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_i\left(t\right))-F_i(t,{\Phi }_2,u_i\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 the solution, there exists a constant $M_1$ such 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) Letting $0<\xi <1$ and $m_k\left(t\right),v_k\left(t\right)\in \Omega$, we obtain $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 $c_1>0$, $\rho >1$, $c_2>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 Ref. [50], we find that there exists an optimal control ${\mu }^{*}=({\mu }_1,\cdots ,{\mu }_n)$, minimizing $J\left(\mu \right)$. The proof is completed.  □

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

$$H=\sum _{k=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 nontrivial 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 that

$$\left\{\begin{array}{cc}{\omega }^{*}=0,\hfill & \frac{\partial H}{\partial u_k\left(t\right)}<0,\hfill \\ 0<{\omega }^{*}<0.9,\hfill & \frac{\partial H}{\partial u_k\left(t\right)}=0,\hfill \\ {\omega }^{*}=0.9,\hfill & \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^{*},\dots ,u_n^{*})$ is given by

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

(25)

Proof. 

We apply the necessary conditions in the Pontryagin maximum principle with Hamiltonian function (21) to derive the adjoint variables 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 conduct some numerical simulations to validate and supplement our theoretical results and to study the impact of the parameters on the dynamics of disease in order to find better control strategies. We first verify and supplement the stability results of system (1) by applying the Ode45 function in Matlab to numerically solve differential equations. It is based on the Runge Kutta method, which can efficiently solve ordinary differential equations and rigid differential equations. Next, we depict the numerical results of the optimal control problem by a combination of forward and backward difference approximation in Matlab to solve the optimal control problem. The numerical simulations include three parts. The first part is to verify the main theoretical results (Theorems 1–4) including the epidemic threshold and global stability, shown in Figure 2 and Figure 3. In the second part, we provide discussion and simulations to show the effectiveness of different immunity and quarantine strategies with degree, as shown in Figure 4 and Figure 5. The numerical simulation in the third part demonstrates the effectiveness of the optimal control method that was proven in Section 5, as shown in Figure 6 and Figure 7. We construct a BA (Barabási–Albert) scale-free network [17], 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 a size of $m_0=3$, and during each time step, a new node with $m=3$ edges was added to the network. We note that $\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 2a shows the relationship between $I\left(t\right)$ and infectious rate $\lambda$ on several different networks. Based on the basic reproduction number $R_0$, we find that 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. Figure 2b and

Table 2. Final infection scale $I\left(t\right)$ under different $R_0$ in Figure 2b.

Basic reproduction number $R_0$	0.1734	0.6590	0.9538	5.3671	7.1561	8.9451
Final infection scale $I\left(t\right)$	0	0	0	0.0449	0.0597	0.0724

verify the results of Theorems 3 and 4, i.e., if $R_0<1$, the disease-free equilibrium is globally stable, and if $R_0>1$, then the endemic equilibrium persists and is globally stable, 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. Figure 2b shows that the disease-free equilibrium cannot undergo a Hopf bifurcation.

Figure 3a shows a bifurcation diagram of model system (3). It clearly appears that model system (3) exhibits a forward bifurcation, that is, the disease-free equilibrium $DFE$ is stable if $R_0\le 1$, while if $R_0>1$, the disease-free equilibrium $DFE$ is unstable, and there exists a unique endemic equilibrium $EE$ which is stable. This result is consistent with the result shown in Figure 2b.

We conduct a parameter sensitivity analysis on the basic reproduction number $R_0$, where

$$R_0=\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}.$$

Then, let $\lambda \left(i\right)=\lambda i$, we can obtain

$$\frac{\partial R_0}{\partial \gamma }=-\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)}^2},\frac{\partial R_0}{\partial \lambda }=\frac{1}{\langle k\rangle }\sum _{i=1}^n\frac{\phi \left(i\right)p\left(i\right)\ast i(1-{\delta }_i)N_i^{*}}{{\beta }_i+\gamma +d}.$$

$$\frac{\partial R_0}{\partial {\beta }_i}=-\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)}^2},\frac{\partial R_0}{\partial {\delta }_i}=-\frac{1}{\langle k\rangle }\sum _{i=1}^n\frac{\phi \left(i\right)p\left(i\right)\lambda \left(i\right)N_i^{*}}{{\beta }_i+\gamma +d},i=1,2,\dots ,n.$$

This indicates that the basic reproduction number $R_0$ increases with the decreases in quarantine rate ${\beta }_k$, immunization rate ${\delta }_k$, and recovery rate $\gamma$, and $R_0$ increases with the increase in infection rate $\lambda$.

Without losing generality, let ${\beta }_k\triangleq \beta$, ${\delta }_k\triangleq \delta$. We depict the relationship between the $R_0$ and parameters in Figure 3b–d, which are consistent with the analysis results. Figure 3b shows that $R_0$ is positively correlated with the birth rate b and negatively correlated with the death rate d.

(II) Secondly, we depict the effectiveness of the immune and quarantining control strategies in Figure 4 and Figure 5. In this part, we also discuss the SIQS model (1) on the network with various immunization (proportional immunization strategy $\tilde{\delta }$, target immunization strategy ${\delta }_k^t$, acquaintance immunization strategy ${\delta }_k^{ac}$) schemes based on [51] and define quarantine schemes of infected individuals (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$. In randomly selecting one infected node for isolation, we find the quarantine rate is independent of the degree of the 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 [51]. Introduce an upper threshold $\kappa$, such that all infected nodes with connectivity $k>\kappa$ are priority quarantined, i.e., we define 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 an 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}}$ denotes the average quarantine rate of acquaintance quarantine, where $\widetilde{{\beta }_k^{ac}}={\sum }_k{\beta }_k^{ac}p\left(k\right)$.

Figure 4a shows the $I\left(t\right)$ under different quarantining control strategies; the target quarantine strategy is the most effective in controlling disease. Figure 4b compares the effectiveness of different quarantine strategies. From Figure 4a,b and

Table 3. Final infection scale and ${\lambda }_c$ under different control strategies in Figure 4a,b.

Quarantining Control Strategies	Final Infection Scale $\mathit{I}\mathbf{\left(}\mathit{t}\mathbf{\right)}$	Epidemic Threshold ${\mathit{\lambda }}_{\mathit{c}}$
No quarantine	0.0713	0.0468
Acquaintance quarantine	0.0568	0.0724
Proportional quarantine	0.0424	0.0798
Target quarantine	0.0182	0.1483

, we find that the target quarantine strategy has better effectiveness than others. Figure 4c shows the effectiveness of different immunity strategies; it shows that all three immunization schemes are effective compared to the case without any immunization, and the targeted immunization scheme is more efficient than the proportional scheme discussed.

Figure 5a depicts the average infectious density with respect to quarantine rate under different immunity strategies. It also shows that the target strategy is better than the proportional and acquaintance strategies. We further show the average infectious density with respect to infected rate for the target and acquaintance immunization schemes and quarantine strategies in Figure 5b. We find 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, as shown in Figure 6 and Figure 7. In Figure 6, we depict $I\left(t\right)$ under different control strategies. In order to provide a clearer representation of the findings in Figure 6a, we additionally compute the objective function $J\left(\mu \right)$ values for diverse control approaches, as shown in

Table 4. The values of $J\left(\mu \right)$ under different control strategies.

Optimal Control	Max Control	No Control
$J\approx 1.7817\times {10}^3$	$J\approx 6.5164\times {10}^4$	$J\approx 5.1127\times {10}^3$

. 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, finding 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 6b, we compare, the effects of optimal control, without control, and acquaintance immunity, combining these with multiple quarantine strategies for infectious disease control. Compared with other strategies, the optimal control strategy has the best control effect.

Figure 7a 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 7b 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. It indicates a higher level of infection, resulting in a higher value of ${\mu }_k$.

7. Conclusions

We proposed an SIQS epidemic model with heterogeneous immune and quarantine measures on heterogeneous networks. We performed a detailed mathematical analysis of our system, which showed 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 of system (1) to tackle the optimal quarantine control issue, confirming the presence of optimal solutions. We hope this result can provide some ideas for solving the optimal solution of asymmetric quarantine strategies on heterogeneous networks. We defined different forms of quarantine measures, which are asymmetric quarantine strategies related in degree, 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 a 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 has a significant cost and workload. Therefore, in order to achieve the control goal and reduce the control cost at the same time, the optimal control strategy considering the asymmetry of control measures related in degree and time is shown to be superior to the symmetric control strategies.

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. This paper can contribute to an understanding of the dynamics of infectious disease with 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 disease data will be meaningful work.