Global dynamics of an endemic disease model with vaccination: Analysis of the asymptomatic and symptomatic groups in complex networks

1.   Introduction

It is well known that many infectious diseases have an incubation period. For example, the average incubation period of HIV is 8–9 years. SARS-CoV-2, which causes COVID-19, also has an incubation period, and people infected by the disease appear as asymptomatic or symptomatic and are able to infect susceptible individuals. To reduce the spread of COVID-19, vaccination is a commonly used control measure. Thus, investigating the impact of vaccines on the dynamics of endemic diseases is crucial.

Many endemic compartmental models for COVID-19 have been proposed to study the transmission dynamics since its outbreak. For example, Ma et al. ^[1] proposed an SEIR-type epidemic model that considers the contact distance between the susceptible individuals and the asymptomatic or symptomatic infected individuals. They analyzed the stability of each equilibrium and obtained the threshold values for population influx and contact distance. Simulation of the SARS-CoV-2 pandemic in Germany has been performed with ordinary differential equations in MATLAB, as presented in ^[2]. It revealed vaccination to be an effective means in the fight against the pandemic. Yuan et al. ^[3] discussed an $\rm SEII_{a}HR$ model and obtained that the disease-free equilibrium (DFE) is globally asymptotically stable if $R_0 < 1$, and that the endemic equilibrium (EE) is uniformly persistent if $R_0 > 1$. Chen et al. ^[4] considered both the SEIR model and the SPEIQRD model incorporating two control actions, i.e., vaccination and quarantine, to estimate the spread of the virus and control the number of infected and dead people via controller tuning. Rakshit et al. ^[5] proposed an SIR model with the asymptomatic and symptomatic groups. They simulated case data from the UK, USA and India with their model. Guo et al. ^[6] proposed an SEIAQR model for COVID-19. By using the limit system of the model and Lyapunov function method, they proved that COVID-19 would die out if the basic reproduction number $R_c \le 1$, and that the disease is persistent if $R_c > 1$. Khoshnaw et al. ^[7] discussed an SIUWR model for COVID-19, and some key computational simulations and sensitivity analysis were investigated. Aguilar et al. ^[8] considered a mathematical model that takes asymptomatic carriers into account. Peirlinck et al. ^[9] modeled the epidemiology of COVID-19 by using an SEIIR model. Simulations indicated that encouraging increased hygiene, mandating the use of face masks, restricting travel and maintaining distance could be the most successful strategies to manage the impact of COVID-19 until vaccination and treatment could become available.

Since the mutation of single-stranded RNA viruses is relatively strong, such as those for HIV, HCV, SARS, MERS-CoV, SARS-CoV-2 and so on, multi-strain models have also been intensively proposed to study epidemic dynamics; for examples, see ^{[10,11,12,13,14]} and the references therein. Many authors have discussed the stability of equilibria in their studies, which may help us determine the long-term behaviors of disease dynamics.

The studies mentioned above mainly made homogeneous contact assumptions, i.e., that each individual in the population has the same probability of contact with an infected individual. However, the contact among people is heterogeneous in the real world. For this reason, the notion of complex networks is incorporated. Many endemic models based on networks have been proposed. Huang et al. ^[15] proposed a network-based SIQRS epidemic model. Sun et al. ^[16] constructed an SIS model on networks. They studied the dynamics of the models, including the basic reproduction number and the global asymptotic stability of the DFE and EE. Meng et al. ^[17] established an SEIRV model on heterogeneous networks. In addition, they used the epidemic parameters in Wuhan, China for numerical analysis. Lv et al. ^[18] proposed an SIVS epidemic model based on scale-free networks. Li and Yousef ^[19] proposed a network-based SIR epidemic model with a saturated treatment function. And, they obtained the threshold values $R_0$ and $\hat{R}_0$. The global stability of equilibria strongly rely on both of them. Besides, forward or backward bifurcation will occur at $R_0 = 1$. Yao and Zhang ^[20] developed a two-strain SIS model based on heterogeneous networks. They derived five different threshold values, i.e., $R_0$, $R_1$, $R_2$, $R_{12}$ and $R_{21}$, which are closely related to the stability of equilibria. Yang and Li ^[21] studied a two-strain SIS epidemic model based on complex networks. The DFE $E_0$ is globally asymptotically stable if the basic reproduction number $R_0 < 1$. Otherwise, there exists either a unique strain 1-only or a strain 2-only boundary equilibrium. Cheng et al. ^[22] talked about a network-based SIQS infectious disease model with a nonmonotonic incidence rate. Moreover, endemic models with time delay in networks have also been established ^[23,24].

As the virus mutates, individuals will become infected and remain asymptomatic. However, some will become seriously ill and require medical attention. More and more people will get vaccinated to prevent infection. Furthermore, the spread of the disease begins exhibiting more heterogeneity. To consider the impact of these features on the dynamics of disease transmission, motivated by the above discussions, we propose an SEIR-type endemic model that incorporates direct immunity by vaccination, as well as the conversion from the asymptomatic to the symptomatic group in complex networks. Conditions for the extinction and permanence of the disease are investigated in this paper.

The rest of the paper is organized in the following way. In Section 2, we introduce the network model. In Section 3, we investigate the dynamics of the proposed model. Section 4 is devoted to numerical simulations to verify our theoretical results. Finally, we give the discussion in Section 5.

2.   Models

As reported by the World Health Organization, COVID-19 affects different people in different ways. Some infected people will be asymptomatic and recover without hospitalization, whereas some will experience fever, coughing, difficulty breathing, shortness of breath, chest pain and other serious symptoms. Both asymptomatic and symptomatic infected individuals are infectious. On average, it takes 5–6 days from when someone is infected with the virus for symptoms to show. According to these COVID-19 transmission characteristics, the total populations are divided into five compartments, namely, the susceptible subclass ($S$), the exposed subclass ($E$), the asymptomatic infected subclass ($I_a$), the symptomatic infected subclass ($I$) and the removed subclass ($R$). The transmission process for each individual is shown in Figure 1. Moreover, we make the following assumptions:

$1)$ All infected people have an incubation period. The incidence functions for the asymptomatic and the symptomatic infected people are $\alpha_1 S I_a$ and $\alpha_2 S I$, respectively.

$2)$ The infected people have immunity to COVID-19 after rehabilitation and vaccination.

$3)$ The asymptomatic infected people may become symptomatic or recover from the disease, which is dependent on their physical fitness. We assume that the rate of transition from $I_a$ to $I$ is proportional to their density with the scaling factor $\theta$.

Based on the above assumptions, we propose the following ordinary differential equation endemic model of COVID-19, which incorporates direct immunity by vaccination and the shift from the asymptomatic to symptomatic subclass:

System (2.1) is an SEIR-type endemic model with two different manifestations of infection in the population. The parameter $A$ is the recruitment rate of susceptible individuals, $\mu$ is the natural mortality rate of all individuals, $\nu$ is the vaccination rate of susceptible individuals, $\alpha_1$ denotes the rate of transmission of susceptible individuals to exposed individuals, as induced by asymptomatic individuals, and $\alpha_2$ denotes the rate of transmission of susceptible individuals to exposed individuals, as induced by symptomatic individuals. Exposed individuals can become asymptomatic at a rate of $\eta_1$, or symptomatic at a rate of $\eta_2$ due to their constitutions. Asymptomatic infected individuals will become symptomatic infected individuals at a rate of $\theta$. $\mu_1$ and $\mu_2$ are the fatality rates of the asymptomatic and the symptomatic individuals, respectively. $\gamma_1$ and $\gamma_2$ are the recovery rates of asymptomatic individuals and symptomatic individuals, respectively. To study the effects of contact heterogeneity on endemic dynamics, we view the entire population as a social network, and each individual corresponds to a network node. Moreover, an edge connecting two individuals denotes a potential contact between both. All nodes of the network are assigned to different groups according to their degrees. Specifically, the $k$th group has a degree $k$ for $k\in\mathbb{N}_n \triangleq \{1, 2, ..., n\}$, and $n$ is the maximal degree. Let $N_k$ be the number of all individuals that have a degree $k$. Then, the total number of individuals is $N = N_1+N_2+...+N_n$. To incorporate the property of contact heterogeneity, we consider the network characterized by the degree distribution $P(k)$, which is defined as a randomly selected node with a degree $k$, i.e., $P(k) = N_k/N$. Based on system (2.1), we construct the following network model with $k$ subsystems for $k \in \mathbb{N}_n$:

where $S_k$ represents the number of susceptible individuals with degree $k$, and $E_k$ represents the number of exposed individuals with degree $k$. $J_k$ and $I_k$ are the numbers of asymptomatic and symptomatic individuals with degree $k$, respectively. $R_k$ is the number of recovered individuals with degree $k$. The parameters for systems (2.1) and (2.2) have the same meanings. Nevertheless, there are slightly different units after dimensionless transformation (see Table 1). And, the transmission diagram for system (2.2) is shown in Figure 2.

All subsystems of dynamics in system (2.2) are coupled by functions $\varTheta_1(t)$ and $\varTheta_2(t)$, which are respectively defined by

where $P(j|k)$ represents the conditional possibility of a node with degree $k$ connecting a node with degree $j$. The possibility of a link connecting to a node with degree $j$ is proportional to $jP(j)$. Hence, $P(j|k) = jP(j)/\langle k \rangle$, where $\langle k \rangle = \sum_{j = 1}^n jP(j)$ is the mean degree of the network ^[25]. Then,

In this paper, we study the endemic dynamics with the form of relations given by (2.3) for the functions $\varTheta_1$ and $\varTheta_2$ in uncorrelated networks.

Since the first four equations in system (2.2) are independent of the variable state $R_k(t)$, we can simplify the model as follows:

3.   Dynamics

Lemma 3.1. Let $k \in \mathbb{N}_n$. For any initial data point $S_k(0)$, $E_k(0)$, $J_k(0)$, $I_k(0) > 0$, all of the solutions ($S_k(t)$, $E_k(t)$, $J_k(t)$, $I_k(t)$) of system (2.4) are nonnegative for $t > 0$.

Proof. Assume, on the contrary, that at least one of the subclass values, i.e., $S_k(t)$, $E_k(t)$, $J_k(t)$, $I_k(t)$, is negative with the initial conditions $S_k(0)$, $E_k(0)$, $J_k(0)$, $I_k(0) > 0$ for $k \in \mathbb{N}_n$. According to the continuity, there exists a sufficiently small $\varepsilon > 0$ such that $S_k(t)$, $E_k(t)$, $J_k(t)$, $I_k(t)\ge 0$ for $t \in (0, \varepsilon)$ and $k \in \mathbb{N}_n$. Furthermore, there exist $j \in \mathbb{N}_n$ and the initial time $t_3 \ge \varepsilon > 0$ such that $J_j(t_3) = 0$ and $J_j^{\prime}(t_3) < 0$ or $I_j(t_3) = 0$ and $I_j^{\prime}(t_3) < 0$. Otherwise, $J_k(t), I_k(t) \ge 0$ for $t > 0$ and $k \in \mathbb{N}_n$, which means that there is no such $t_3$. Similarly, there exist $j \in \mathbb{N}_n$ and the initial time points $t_1$ and $t_2$ such that $S_j(t_1) = 0$, $S_j^{\prime}(t_1) \le 0$ and $E_j(t_2) \le 0$, $E_j^{\prime}(t_2) = 0$ for $t_1, t_2 \ge \varepsilon > 0$.

Case 1: $t_1\le t_2, t_3$. Here $t_2$ and $t_3$ may not exist, but $t_3$ must exist. We have that $E_j(t_1), J_j(t_1), I_j(t_1), \varTheta_1(t_1), \varTheta_2(t_1) \ge 0$. Substituting $t_1$ into the first equation of system (2.4) yields $S_k'(t_1) = A > 0$. Hence, this contradicts the hypothesis.

Case 2: $t_2 \le t_1, t_3$. Here $t_1$ and $t_3$ may not exist, but $t_2$ must exist. We have that $S_j(t_2), J_j(t_2), I_j(t_2), \varTheta_1(t_2), \varTheta_2(t_2) \ge 0$. Substituting $t_2$ into the second equation of system (2.4) yields $E_j'(t_2) = \alpha_1 j S_j(t_2) \Theta_1(t_2)+\alpha_2 j S_j(t_2) \Theta_2(t_2) \ge 0$. Hence, this contradicts the hypothesis.

Case 3: $t_3 \le t_1, t_2$. Here $t_1$ and $t_2$ may not exist, but $t_3$ must exist. If $J_j(t_3) = 0$, we have that $S_j(t_3), E_j(t_3), I_j(t_3), \varTheta_1(t_3), \varTheta_2(t_3) \ge 0$. Substituting $t_3$ into the third equation of system (2.4) yields $J_j'(t_3) = \eta_1E_j(t_3) \ge 0$. Hence, this contradicts the hypothesis. If $I_j(t_3) = 0$, we have that $S_j(t_3), E_j(t_3), J_j(t_3), \varTheta_1(t_3), \varTheta_2(t_3) \ge 0$. Substituting $t_3$ into the last equation of system (2.4) yields $I'(t_3) = \eta_2E_j(t_3)+ \theta J_j(t_3) \ge 0$. Hence, this contradicts the hypothesis.

Similarly, we can prove that $S_k(t)\geq0$. This completes the proof. □

Lemma 3.2. The solutions of (2.4) satisfy that $0\leq S_{k}\leq \frac{A}{\mu}$, $0\leq E_{k}\leq \frac{A}{\mu}$, $0\leq J_{k}\leq \frac{A}{\mu}$ and $0\leq I_{k}\leq \frac{A}{\mu}$ for all $k \in \mathbb{N}_n$.

Proof. Let $N_k = S_{k}+ E_{k}+ J_{k}+ I_{k}$ for $k\in\mathbb{N}_{n}$. Adding the equations in system (2.4), we obtain

By integration, it follows that

where $N_k(0)$ represents the total population with degree $k$ at time $t = 0$. Therefore,

This, together with Lemma 3.1, readily imply the lemma. □

System (2.4) always has the DFE

where

Assume that there exists the EE

where

$k \in \mathbb{N}_n$. Solving system (3.1), we obtain that $S_k^*, E_k^*, J_k^*$ and $I_k^*$ are functions with respect to $\varTheta_1$ and $\varTheta_2$ for $k \in \mathbb{N}_n$, namely,

Substituting (3.3) into (2.3) gives

From (3.4) and (3.5), we observe

Noting (3.4)–(3.6), we define

Then, (3.4) and (3.5) are respectively equivalent to

Clearly, $\varTheta_1 = 0$ and $\varTheta_2 = 0$ are solutions for (3.9), which correspond to the DFE $P_0$. Assume that there exist nontrivial solutions $0 < \varTheta_1 < 1$ and $0 < \varTheta_2 < 1$ for (3.9). Then, we consider derivatives for the functions (3.7) and (3.8). Hence, we have

for $\varTheta_1 > 0$ and $\varTheta_2 > 0$. Owing to $f_1(0) = 0$ and $f_2(0) = 0$, it is guaranteed that the equations given by (3.9) have nontrivial solutions between 0 and 1 if

Combining (3.10) with (3.7) and (3.8) implies that

and

Define the basic reproduction number of system (2.4) as

Hence, there exists a unique positive EE $P^*$ with $0 < \varTheta_1 < 1$ and $0 < \varTheta_2 < 1$ for system (2.4) if

where

Remark 3.1. There always exists a DFE $P_0$ for system (2.4). Besides, given that $\lim\limits_{\varTheta_1 \to +\infty} f_1(\varTheta_1) = \lim\limits_{\varTheta_2 \to +\infty} f_2(\varTheta_2) = -\infty$, and considering the continuity of functions $f_1(\varTheta_1)$ and $f_2(\varTheta_2)$, there exist solutions $\varTheta_1 > 0$ and $\varTheta_2 > 0$ for (3.9) if and only if $R_0 > 1$. Namely, system (2.4) has a unique positive EE $P^*$ if and only if $R_0 > 1$.

Theorem 3.1. The DFE $P_0$ for system (2.4) is locally asymptotically stable if $R_0 < 1$, and unstable if $R_0 > 1$.

Proof. The Jacobian matrix for system (2.4) is

where $a = \mu+\nu$, $b = \mu+\eta_1+\eta_2$, $c = \mu+\mu_1+\gamma_1+\theta$, $d = \mu+\mu_2+\gamma_2$, $k \in \mathbb{N}_n$. Obviously, there are $n$ eigenvalues equaling to $-\mu-\nu$ for the matrix (3.11), and the remaining eigenvalues are determined by

$k \in \mathbb{N}_n$. For (3.12), we add the column $2n+1$, multiplied by $-\frac{\alpha_1}{\alpha_2}$, to column $n+1$; column $2n+2$, multiplied by $-\frac{\alpha_1}{\alpha_2}$, to column $n+2$ $\cdots$; and column $3n$, multiplied by $-\frac{\alpha_1}{\alpha_2}$, to column $2n$. Hence, (3.12) becomes

$k \in \mathbb{N}_n$. For (3.13), we add the row $n+1$, multiplied by $\frac{\alpha_1 d+\alpha_2 \theta}{\alpha_2 c}$, to row $2n+1$; row $n+2$, multiplied by $\frac{\alpha_1 d+\alpha_2 \theta}{\alpha_2 c}$, to row $2n+2$; $\cdots$ and row $2n$, multiplied by $\frac{\alpha_1 d+\alpha_2 \theta}{\alpha_2 c}$, to row $3n$. Then, (3.13) is converted to

$k \in \mathbb{N}_n$. For (3.14), we add the column $n+1$, multiplied by $\frac{\eta_1}{c}$, to column $1$; column $n+2$, multiplied by $\frac{\eta_1}{c}$, to column $2$; $\cdots$ and column $2n$, multiplied by $\frac{\eta_1}{c}$, to column $n$. Then, (3.14) is converted to

$k \in \mathbb{N}_n$. Obviously, there are $n$ eigenvalues equaling to $-\mu-\nu_1-\gamma_1-\theta$ for (3.15). The remaining eigenvalues are determined by

$k \in \mathbb{N}_n$. For (3.16), we add the column $n+2$, multiplied by $-\frac{P(1)}{2P(2)}$, to column $n+1$; column $n+3$, multiplied by $-\frac{2P(2)}{3P(3)}$, to column $n+2$; $\cdots$ and column $2n$, multiplied by $-\frac{(n-1)P(n-1)}{nP(n)}$, to column $2n-1$. Then, (3.16) is converted to

$k \in \mathbb{N}_n$. For (3.17), we add the row $1$, multiplied by $\frac{(\alpha_1 d+\alpha_2 \theta)\eta_1+\alpha_2 c \eta_2} {\alpha_2 b c}$, to row $n+1$; row $2$, multiplied by $\frac{(\alpha_1 d+\alpha_2 \theta)\eta_1+\alpha_2 c \eta_2} {\alpha_2 b c}$, to row $n+2$; $\cdots$ and row $n$, multiplied by $\frac{(\alpha_1 d+\alpha_2 \theta)\eta_1+\alpha_2 c \eta_2} {\alpha_2 b c}$, to row $2n$. Then, (3.17) is converted to

$k \in \mathbb{N}_n$. Obviously, there are $n$ eigenvalues equaling to $-\mu-\eta_1-\eta_2$ for (3.18); the remaining eigenvalues are determined by

$k \in \mathbb{N}_n$. For (3.19), we add the row $1$, multiplied by $\frac{P(1)}{2P(2)}$, to row $2$; row $2$, multiplied by $\frac{2P(2)}{3P(3)}$, to row $3$; $\cdots$ and row $n-1$, multiplied by $\frac{(n-1)P(n-1)}{nP(n)}$, to row $n$. Hence, (3.19) is converted to

$k \in \mathbb{N}_n$. Clearly, there are $n-1$ eigenvalues equaling to $-\mu-\mu_2-\gamma_2$ for (3.20), and the last eigenvalue is given by

Hence, all eigenvalues are negative if $R_0 < 1$, and there exists a positive eigenvalue if $R_0 > 1$. Therefore, the DFE $P_0$ for system (2.4) is locally asymptotically stable if $R_0 < 1$, and unstable if $R_0 > 1$.

Remark 3.2. Define

Then, from Theorem 3.3, we can derive the threshold-type dynamics of the DFE of system (2.4) with respect to $A^*$, namely, the DFE is locally asymptotically stable if $A < A^*$, and unstable if $A > A^*$.

Theorem 3.2. Let

If $A \le A^g$, then the DFE $P_0$ for system (2.4) is globally asymptotically stable.

Proof. Let $h(x) = x-1-\ln{x}$. It is easy to verify that the function $h(x)$ is nonnegative. We define the Lyapunov function

Obviously, $V \ge 0$, and the equality holds if and only if $E_k = J_k = I_k = 0$ and $S_k = \frac{A}{\mu+\nu}$ for $k \in \mathbb{N}_n$. Differentiating $V$ with respect to $t$, we derive

Then, $V^{\prime}(t) \le 0$ is guaranteed if $A \le A^g$, and the equality holds if and only if $E_k = J_k = I_k = 0$ and $S_k = \frac{A} {\mu+\nu}$, $k \in \mathbb{N}_n$. Hence, by LaSalle's invariance principle ^[26], $P_0$ is globally asymptotically stable. This completes the proof. □

Remark 3.3. Note that $A < A^*\Leftrightarrow R_0 < 1$. A direct computation yields $A^* > A^g$. Then, we have

Hence, condition $A \le A^g$ is stronger than condition $R_0 < 1$.

For the global asymptotic stability of EE, we assume that

Consequently, system (2.4) becomes

where $\alpha = \alpha_1 = \alpha_2$, $\eta = \eta_1+\eta_2$, $\varTheta = \varTheta_1+\varTheta_2$, $L_k = I_k+J_k$. And, the invariant set for system (3.22) is

The EE of system (3.22) is $\widetilde{P}^* = (S_k^*, E_k^*, L_k^*)$, $k \in \mathbb{N}_n$, where

Owing to (3.21), we obtain the basic reproduction number

Theorem 3.3. The EE $\widetilde{P}^*$ of system (3.22) is globally asymptotically stable in $\widetilde{\Omega}$ if $R_0 > 1$ and the matrix $\left(\frac{\alpha k k' P(k')}{\langle k \rangle}\right)_ {n \times n}, k = 1, 2, \cdots, n$, $k' = 1, 2, \cdots, n$ is irreducible.

Proof. Define

Hence, $V_k \ge 0$, and the equality holds if and only if $S_k = S_k^*$, $E_k = E_k^*$ and $L_k = L_k^*$. Differentiating $V_k$, we obtain

Let $a_{kk^{\prime}} = \alpha k S_k^* \sum_{k^\prime = 1}^{n}\frac{k^\prime P(k^\prime) L_{k^\prime}^*}{\langle k \rangle}$, $G_k(L_k) = -\frac{L_k}{L_k^*}+\ln{\frac{L_k}{L_k^*}}$ and $H(x) = -h(x) = 1+\ln{x}-x$. We have

Therefore, $a_{kk'}$, $G_k$, $F_{kk'}$ and $V_k$ satisfy the assumptions of Theorem 3.1 and Corollary 3.3 in Ref. ^[27]. For the Lyapunov function $V^* = \sum_{k = 1}^{n} c_k V_k$, we have that $\frac{dV^*}{dt} \le 0$ for $(S_1, E_1, L_1, S_2, E_2, L_2, \cdots, S_n, E_n, L_n) \in \widetilde{\Omega}$. We can show that the only compact invariant set $\frac{dV^*}{dt} = 0$ is a singleton $\widetilde{P}^*$. By LaSalle's invariance principle, we conclude that $\widetilde{P}^*$ is globally asymptotically stable in $\widetilde{\Omega}$ if $R_0 > 1$. This completes the proof.

Remark 3.4. According to Theorem 3.5, the EE will be persistent if $R_0 > 1$ and the matrix $\left(\frac{\alpha k k' P(k')}{\langle k \rangle}\right)_{n \times n}$ is irreducible for $k$, $k^\prime \in \mathbb{N}_n$. To compare the threshold-type dynamics with respect to the reproduction number, we further need the irreducible assumption. This restriction is just a technical requirement in computation.

4.   Simulations

In this section, we present numerical simulations to verify our theoretical results. We set the initial values $S_k(0) \in (0.5, \, 0.6)$, $E_k(0) \in (0, \, 0.1)$, $J_k(0) \in (0, \, 0.1)$ and $I_k(0) \in (0, \, 0.1)$ for $k \in \mathbb{N}_n$ in all simulations. We consider a scale-free network with a degree distribution of $P(k) = (\gamma-1) m^{\gamma-1} k^{-\gamma}$, where $m$ is the smallest degree of node in a scale-free network; $\gamma$ stands for the power law exponent. Here, we applied $m = 2, \gamma = 3$ and $\mathbb{N}_n = \{1, 2, \cdots, 100\}$. Consequently, we obtained that $\langle k \rangle = 6.5399$ and $\langle k^2 \rangle = 20.7495$.

In Figure 3, the parameters were selected as $A = 0.02$, $\mu_1 = 0.01$, $\mu_2 = 0.02$, $\eta_1 = 0.4$, $\eta_2 = 0.2$, $\gamma_1 = 0.4$, $\gamma_2 = 0.6$, $\nu = 0.01$, $\alpha_1 = 0.02$, $\alpha_2 = 0.04$, $\theta = 0.02$ and $\mu = 0.01$. Subsequently, we have that $R_0 = 0.1666 < 1$, and it follows from Theorem 3.1, i.e., that the DFE is locally asymptotically stable. Moreover, according to Theorem 3.2, we have that $A^g = 0.0993$, which implies that $A < A^g$. This means that the DFE is globally asymptotically stable according to our theoretical analysis. From Figure 3, it can be seen that $S_k(t) \to S_k^0(t) = 1$, and $E_k(t)$, $J_k(t)$, $I_k(t)$, $\varTheta_1(t)$, $\varTheta_2(t)\to 0$ as $t\to\infty$, where $k = 20, 40, 60, 80, 100$, respectively. This supports the stability of the DFE.

In Figure 4, the parameters were selected as $A = 0.02$, $\mu_1 = \mu_2 = 0.01$, $\eta_1 = 0.4$, $\eta_2 = 0.2$, $\gamma_1 = \gamma_2 = 0.4$, $\nu = 0.01$, $\alpha_1 = \alpha_2 = 0.2$, $\theta = 0.02$ and $\mu = 0.01$. Subsequently, we have that $R_0 = 1.4861 > 1$, and it follows from Theorem 3.3 that the endemic equilibrium $\widetilde{P}^*$ of system (3.22) is globally asymptotically stable. Namely, the disease will be endemic. Figure 4 shows that $S_k(t)$, $E_k(t)$, $J_k(t)$, $I_k(t)$, $\varTheta_1(t)$ and $\varTheta_2(t)$ will persist for a sufficiently large $t$, where $k = 20, 40, 60, 80, 100$. This supports the stability of the EE.

In Figures 5 and 6, we consider the different recruitment rates in groups with different degrees. Let $A_k$ be the recruitment rate of the group with degree $k$. The parameters were selected to be the same as those applied in Figures 3, and $A_k$ was randomly selected in $\left(0, A^g\right]$ for $k \in \mathbb{N}_n$, where $A^g = 0.0993$. The simulation results are shown in Figures 5. It can be seen that $S_k(t)$ will persist, and $E_k(t)$, $J_k(t)$, $I_k(t)$, $\varTheta_1(t)$, $\varTheta_2(t) \to 0$ for a sufficiently large $t$ and $k = 20, 40, 60, 80, 100$. Therefore, it seems that, as long as all values of $A_k$ are below the threshold for $k \in \mathbb{N}_n$, the disease will die out.

In Figures 6, the parameters were selected to be the same as those applied in Figures 4 and $A_k$ was randomly chosen in $\left(A^*, A_0\right]$ for $k \in \mathbb{N}_n$, where $A^* = 0.0135$ and $A_0 = 0.6710$. It can be seen that $S_k(t)$, $E_k(t)$, $J_k(t)$, $I_k(t)$, $\varTheta_1(t)$ and $\varTheta_2(t)$ will persist for a sufficiently large $t$, where $k = 20, 40, 60, 80, 100$. Therefore, it seems that, as long as all values of $A_k$ are above the threshold for $k \in \mathbb{N}_n$, the disease will also be endemic.

In Figures 7 and 8, we show the time evolution diagrams for $\Theta_1$ and $\Theta_2$ with respect to $t$ under two different population influx patterns, i.e., the fixed population inflow and the random population influx. The fixed population influx pattern means that the recruitment rates of different groups with different degrees are the same, while the random population influx pattern means that the recruitment rates of different groups are different. We applied random values to depict random recruitment rates in the simulations. Note that the stability of the DFE and EE depends on the strength of the relationship between the recruitment rate and threshold values $A^g$ and $A^*$. Therefore, once the other parameters are fixed, we can obtain the value ranges for $A$ and $A_k, k \in \mathbb{N}_n$ by taking $A^g$ and $A^*$ into account. Figure 7 (left) (Figure 8 (left)) indicates that, if $A \leq A^g$ ($A_k \leq A^g, k \in \mathbb{N}_n$), the disease will be extinct, and Figure 5 (right) (Figure 8 (right)) indicates, that if $A > A^g$ ($A_k > A^*, k \in \mathbb{N}_n$), the disease will be endemic. This is in line with our expectations.

In Figures 9 and 10, we consider the impact of the network structure on the dynamics, and the parameters were fixed as $A = 0.1$, $\mu_1 = \mu_2 = 0.01$, $\eta_1 = 0.4$, $\eta_2 = 0.2$, $\gamma_1 = \gamma_2 = 0.4$, $\nu = 0.01$, $\alpha_1 = \alpha_2 = 0.02$, $\theta = 0.02$ and $\mu = 0.01$. In Figure 9, we consider scale-free networks with the degree distribution $P(k) = (\gamma-1) m^{\gamma-1} k^{-\gamma}$ for $k \in \mathbb{N}_n$, $m = 2$ and $\gamma = 2.1, \ 2.3, \ 2.5, \ 2.7, \ 2.9$, respectively. In Figure 10, we consider Poisson networks with the degree distribution $P(k) = e^{-c}c^k/k!$ for $k \in \mathbb{N}_n$ and $c = 3, \ 6, \ 9, \ 12, \ 15$, respectively. Obviously, we can observe that $\varTheta_1(t)$ and $\varTheta_2(t)$ tend to 0 for a sufficiently large $t$ when $\gamma$ is large in Figure 9. In contrast, $\varTheta_1(t)$ and $\varTheta_2(t)$ approach 0 for a sufficiently large $t$ when $c$ is small in Figure 10. In summary, the network structure does affect the spreading dynamics.

5.   Conclusions

Complex networks, such as completely random networks ^[28], small-world networks ^[29] and scale-free and Poisson networks ^[30], have been frequently used to model the spread of an epidemic disease. Affected by occupation and geographical location, the contact among the population cannot always be modeled as a uniform collision for COVID-19; also, the epidemic disease transmission is usually heterogeneous. In order to investigate the effect of contact heterogeneity, we proposed an endemic mathematical model that incorporates direct immunity by vaccination, as well as the shift from the asymptomatic to the symptomatic subclass, by applying the idea of a compartmental model in a scale-free network. With the help of eigenvalues of the Jacobian matrix, the Lyapunov function method and LaSalle's invariance principle, we evaluated the dynamics of the proposed model. The results were obtained as follows:

$\bullet$ If $R_0 < 1$ or $A < A^*$, the DFE is locally asymptotically stable.

$\bullet$ If $A \le A^g$, the DFE is globally asymptotically stable.

$\bullet$ If $R_0 > 1$ and the matrix $\left(\frac{\alpha k k' P(k')} {\langle k \rangle}\right)_{n \times n}$, $k = 1, 2, \cdots, n, k' = 1, 2, \cdots, n$ is irreducible, then the EE of (2.4) with $\alpha_1 = \alpha_2$, $\mu_1 = \mu_2$ and $\gamma_1 = \gamma_2$ is globally asymptotically stable.

In the past years, ordinary differential equation compartmental models have been widely applied to describe the spread of COVID-19 (see ^{[1,2,3,4,5,6,7,8,9]} for example). Taking contact heterogeneity into account, we have generalized the results of the dynamics for the ordinary differential equations to the networks. We also performed simulations to validate our theoretical results. In the simulations, we mainly show two different scenarios, namely, disease persistence and disease disappearance. So, we selected two different sets of data to simulate. Moreover, we found that different subgroups with different degrees may have different recruitment rates. Again, we selected two more different sets of data to simulate. Based on the result presented in this paper, we know that the recruitment is crucial. Hence, the investigation of the recruitment rates in different groups is an important direction that is significant to better connect with reality.

The paper indicates that the disease may be under control if $R_0 < 1$. So, we can adjust the values of parameters to reduce the value of $R_0$. However, some values of parameters, including the natural mortality rate, transmission rate, recovery rate, disease mortality rate, incubation period duration and the rate of shift from the asymptomatic to the symptomatic group, are inherent attributes of diseases or nature that cannot be changed in the short term. Once these values are determined, we can control the spread of the disease by adjusting the value of the vaccination rate $\nu$ and/or the recruitment rate $A$. From the result of this paper, we see that accelerating vaccination and reducing travel and contact are appropriate strategies for controlling the spread of the disease. Our study can give theoretical perspective that can facilitate understanding of the transmission mechanism for COVID-19, and the result may provide some reasonable suggestions that can be used by the policy-makers to control the disease.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

The authors would like to express their sincere gratitude to the reviewers for providing valuable comments and suggestions that have significantly contributed to the enhancement of the quality of this manuscript. This work was supported by the Youth Talent Promotion Project of Henan Province (No. 2019HYTP035), China Postdoctoral Science Foundation (No. 2018M630824) and Natural Science Foundation of Henan Province (No. 232300420357).

Conflict of interest

The authors declare that there is no conflict of interest.