key: cord-0943270-ysanzw8l
authors: Guan, Gui; Guo, Zhenyuan
title: Stability behavior of a two-susceptibility SHIR epidemic model with time delay in complex networks
date: 2021-08-30
journal: Nonlinear Dyn
DOI: 10.1007/s11071-021-06804-6
sha: 696db0726872887f50324d70aca69757c904da28
doc_id: 943270
cord_uid: ysanzw8l

Taking two susceptible groups into account, we formulate a modified subhealthy-healthy-infected-recovered (SHIR) model with time delay and nonlinear incidence rate in networks with different topologies. Concretely, two dynamical systems are designed in homogeneous and heterogeneous networks by utilizing mean field equations. Based on the next-generation matrix and the existence of a positive equilibrium point, we derive the basic reproduction numbers [Formula: see text] and [Formula: see text] which depend on the model parameters and network structure. In virtue of linearized systems and Lyapunov functions, the local and global stabilities of the disease-free equilibrium points are, respectively, analyzed when [Formula: see text] in homogeneous networks and [Formula: see text] in heterogeneous networks. Besides, we demonstrate that the endemic equilibrium point is locally asymptotically stable in homogeneous networks in the condition of [Formula: see text] . Finally, numerical simulations are performed to conduct sensitivity analysis and confirm theoretical results. Moreover, some conjectures are proposed to complement dynamical behavior of two systems.

and R 2 0 which depend on the model parameters and network structure. In virtue of linearized systems and Lyapunov functions, the local and global stabilities of the disease-free equilibrium points are, respectively, analyzed when R 1 0 < 1 in homogeneous networks and R 2 0 < 1 in heterogeneous networks. Besides, we demonstrate that the endemic equilibrium point is locally asymptotically stable in homogeneous networks in the condition of R 1 0 > 1. Finally, numerical simulations are performed to conduct sensitivity analysis and confirm theoretical results. Moreover, some conjectures are proposed to complement dynamical behavior of two systems.

From historical events, infectious diseases, such as cholera [1] , malaria [2] , influenza [3] and COVID-19 [4] , pose a huge threat to the public health all over the world. Fortunately, mathematical modeling is regarded as a powerful method to characterize the spreading mechanism of the epidemic, which makes it convenient for us to study the dynamics of disease propagation in epidemiology significantly [5] [6] [7] [8] [9] .

According to the classical theory of epidemic dynamics, many scholars progressively apply and extend compartmental epidemic models to research transmission dynamics in a diverse range of fields [10] [11] [12] [13] [14] [15] . In recent years, due to the nonuniformity of spread in a population, the thought of classification for susceptible, infected or other state groups is embedded in the modeling of propagation phenomena in Refs. [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] . Dating back to 1976, a multiple groups model proposed by Lajmanovich and Yorke [29] characterizes the epidemiological features of gonorrhea appropriately. In a general way, the specific group is further divided into several disjoint subgroups in the multi-group model, which enriches epidemiological states of individuals. As a result, the application of classifying thought in defining compartmental states can promote the established propagation model more practical. Throughout previous papers, the classification for a certain group is mostly attributed to genetic variation (e.g., gender). For instance, Hyman and Li [17] , respectively, formu-late thresholds for the spread of infectious disease in the differential susceptibility model, staged progression model, differential infectivity model and so on. Besides, they [18] also study SIR epidemic models with differential susceptibility by classifying the susceptible into various subgroups. Additionally, Wang et al. [20, 21] divide the infected compartment into two sub-compartments and discuss the spreading dynamics of a sexually transmitted disease model when lowrisk and high-risk infected individuals coexist. Considering the heterogeneity of host population, Jin et al. [23] investigate the disease transmission dynamics by establishing a general multi-group epidemic model. In view of the coinfection of two strains, Ruan et al. [25] calculate the basic reproduction number and study the threshold dynamics of a diffusive SIS epidemic model where infected individuals are split into the infected with strain one and strain two based on Ref. [26] .

In Ref. [30] , an SIR epidemic model with consideration of birth, death and two susceptibility is proposed and studied in heterogeneous networks as follows:

(1.1)

The total population is divided into four disjoint compartments including the 1st susceptible group, the 2nd susceptible group, infected individuals and recovered individuals. Their densities are, respectively, denoted by S 1k (t), S 2k (t), I k (t) and R k (t) with degree k (k = 1, 2, . . . , n) at time t. Besides, let T k (t) = S 1k (t) + S 2k (t) + I k (t) + R k (t) be the overall density of individuals with degree k at time t. Here, (1 − T k (t)) represents the density of empty nodes which can generate newborns belonging to the 1st or 2nd susceptible group at the certain rate b i (i = 1, 2). μ is a natural death probability which is identical for all individuals. β 1 or β 2 is the transmission coefficient at which a 1st or 2nd susceptible individual is infected with this disease by getting in contact with infected individuals. Meanwhile, infected individuals can turn into the recovered state at a recovery rate γ . Besides, Θ is known as the probability with which any chosen edge is linked to an infected individual. To be specific, in the uncorrelated networks, Θ can be written as

where p(k ) is the degree distribution and k = k kp(k)

is the average degree of the network. For system (1.1), Yuan et al. [30] obtain the basic reproduction number, analyze the stability of two equilibria and give the effectiveness of control strategies. Nevertheless, the epidemiologically meaningful time delay is ignored in the spreading process of contagious disease in their modeling, which deviates the realistic situation. Strictly speaking, hysteresis involved in virus production indeed exists because some time is needed for the maturity of the virion after the cells of individuals have catched the virus [31] . That is to say, the infected individuals become infectious and further transmit the infection after a certain period of time. Based on this situation, it's of biological significance to bring in the time delay in the transition process due to the existence of latent infection. Recently, lots of researchers are passionate about adding the time delay into the epidemic model which consists of a coupled system of delay differential equations in Refs. [32] [33] [34] [35] [36] [37] [38] . Actually, the introduction of time delay may lead to the change in stability of equilibrium points of a dynamical system, which puts forward a challenge to the dynamical analysis of the delayed model. For example, the existence of Hopf bifurcations at various equilibria in a delayed predator-prey model is explored by Xu in Ref. [32] . Moreover, Zhu et al. [34] study the local and global asymptotic stabilities of equilibrium points of a rumor spreading model with and without time delay, respectively. Considering both avian population and human population, Kang et al. [36] employ two discrete time delays τ 1 and τ 2 to delineate the delayed process in state transitions. In consequence, it is highly meaningful for us to incorporate the time delay into the epidemic model (1.1) additionally.

Similar to Ref. [30] , most of the existing results focusing on epidemic models usually suppose that the interaction term between the susceptible and infected individuals satisfies the bilinear form according to the law of mass action [39] . Once the bilinear incidence rate is used in the modeling, the number of patients linearly increases in the infection process, which isn't appropriate in the situation of huge numbers of infected individuals [40] . In fact, if the epidemic is severe enough, the information about the prevalence of disease will impel individuals to take prevention measures to avoid the infection [41] . In consequence, it's a pity that some researchers don't think over the behavioral changes of individuals due to the psychological effect in disease transmission process. In 1978, to characterize saturation phenomena for mass infected individuals, Capasso and Serio [42] generalize the Kermack-McKendrick deterministic model by employing a saturated incidence rate Sg(I ) where g(I ) = β I 1+a I is a nonlinear bounded function. As explained in Refs. [43, 44] , the incidence function g(I ) gradually reaches a saturation level with the scale of infected individuals I increasing.

Noticing the insightful effect of nonlinear incidence on epidemic dynamics, lots of authors [40, [45] [46] [47] [48] are devoted to studying the nonlinear dynamics of various epidemic models incorporated with saturated incidence rate. For instance, Zhu et al. [48] perfectly explore the stability of equilibrium points and the effectiveness of control schemes in a delayed SIS epidemic model along with nonlinear incidence rate. Furthermore, Li and Yousef [49] analytically and numerically research the bifurcation behavior of a network-based SIR epidemic model with saturated treatment function which is analogous to the nonlinear type of incidence rate [50] in some sense. To address the mentioned deficit in Ref. [30] , we intend to adopt the nonlinear incidence rate by introducing a psychological factor a, which avoids the unbounded contact rate in the epidemic model.

What's more, the application of complex networks to epidemic modeling gives rise to a wave of research in academia for decades. With respect to the consideration of network topology, more and more scholars are in favor of investigating the dynamics of propagation model in both homogeneous and heterogeneous networks [51] [52] [53] [54] [55] . In Ref. [51] , Xia et al. present an improved SEIR model with hesitating mechanism and analyze the spreading threshold in homogeneous and heterogeneous networks, respectively. As Zhu and Guan [53] present, the complexity of the network structure can result in the difference of spreading threshold of the disease in magnitude. In addition, the basic reproduction number and the rumor-free equilibrium point of I2S2R rumor propagating model with two rumors are explored in homogeneous networks by Zhang and Zhu [55] . And they are also devoted to proving the stability of the trivial equilibrium point, discussing the global attractivity of the positive equilibrium point and investigating the permanence of system in heterogeneous networks. As the society developing greatly, the connections among people are increasingly convenient and frequent. This situation drives us to take advantage of complex networks to capture the features of social network in reality. According to insights from studying complex networks, homogeneous and heterogeneous networks can be selected, respectively, as the underlying network to investigate the influence of network structure on disease transmission.

As far as we know, however, there are few research results on dynamical analysis of the delayed multigroup epidemic model with nonlinear incidence rate and different topological structures of social network. Motivated by foregoing literature, such as [19, 30] , in which Yuan et al. study the stability of the SIR epidemic model with differential susceptibility or infectivity, and [48, 53] , in which Zhu et al. incorporate the saturated incidence rate with time delay into propagation model, we establish a delayed epidemic model along with two susceptible groups and nonlinear incidence rate in complex networks. Emphasizing the influence of network structure on disease propagation, we make efforts to investigate threshold dynamics of a two-susceptibility epidemic model with time delay and nonlinear incidence rate in both homogeneous and heterogeneous networks. What's more, to make recommendations for control measures, we pay attention to the impacts of time delay, nonlinear incidence rate and network structure on the transmission of infectious disease.

To make this epidemic model with two susceptible groups more sensible, assume that personal fitness level results in differences between susceptible individuals. In other words, the susceptibility of individuals to infectious disease depends on the personal fitness level. Taking the health level as the classification criterion of susceptible population, we introduce the subhealthy and healthy compartments to characterize the nonhomogeneous structure for susceptible individuals. Hence, states of individuals cover the subhealthy (S), the healthy (H ), the infected (I ) and the recovered (R). Based on the practical situation, it's further supposed that the subhealthy are more likely to catch the infection than the healthy herein. As for the modeling of disease propagation, the network topology and state-transition rules of nodes need to be considered primarily. Since the choice of the network topological structure is dis-cussed above, the expression of interaction rules in the spreading process is briefly presented below.

(1) In virtue of empty nodes, newborns with two levels of health enter the network and all individuals naturally emigrate the network owing to the death.

(2) Upon contacting infectious vectors that carry the pathogen, a subhealthy or healthy individual will turn into an infected individual with a certain probability. Importantly, a nonlinear incidence rate can reflect the crowding effect of infected individuals and inhibitory measures taken by susceptible individuals when the diffusion of disease is especially serious [56] .

(3) As a matter of fact, the time delay plays a major role in the epidemic model because the incubation period of disease indeed exists on account of the latency in a vector. Namely, some time is needed for infectious agents developing in the vector, after which infected vectors become infectious and can transmit the infection to humans.

(4) With the aid of modern treatment, infected individuals can be cured and become the recovered state with a certain rate.

Based on the above analysis, we, respectively, establish the delayed two-susceptibility epidemic model with nonlinear incidence rate in homogeneous networks in Sect. 2 and heterogeneous networks in Sect. 3. Basic properties of solutions and the threshold of disease diffusion are analytically derived. Furthermore, we prove the stability of equilibrium points of two dynamical systems in detail. In Sect. 4, quantities of numerical experiments are carried out to verify the correctness of obtained theoretical results. Finally, the paper ends with conclusions and discussions in Sect. 5.

At first, the delayed system about disease diffusion is considered on the topological structure of homogeneous networks where all nodes are regarded as equivalent statistically. Let S(t), H (t), I (t) and R(t) represent the average densities of subhealthy, healthy, infected and recovered nodes at time t, respectively. The mean field equations of the SHIR epidemic model in homogeneous networks are composed of a set of delay differential equations as follows:

denotes the density of newborn susceptible nodes generated by empty nodes with the certain constant rate b i (i = 1, 2). The psychological factor a characterizes the behavioral changes resulted from the crowding effect of infected individuals during a peak period of epidemic situation. Latent period is shown by the average time delay τ of disease propagation in the process of infection from the subhealthy and healthy to infected individuals. All parameters in our epidemic model are assumed to be positive. The initial conditions for system (2.1) are given by S(ϑ), H (ϑ), I (ϑ), R(ϑ) = ϕ 1 (ϑ), ϕ 2 (ϑ), ϕ 3 (ϑ), ϕ 4 (ϑ) which satisfy

As the sum of S(t), H (t), I (t) and R(t), the total population size at time t can be expressed by T (t). Add up all equations of system (2.1) and obtain the following differential equation:

Making allowances for the above equation, we transform system (2.1) into the limit system:

where T * = b 1 +b 2 b 1 +b 2 +μ . Now, it suffices to study system (2.3) detailedly instead of system (2.1) when exploring the long-time behavior for the solutions of our model. Moreover, the fourth equation of system (2.3) is decoupled from the equations for S(t), H (t) and I (t). Hence, it's natural to make the limit system (2.3) be reduced as the following system:

(2.4)

For the sake of simplicity, it's high time to study system (2.4) in place of original system (2.1) sufficiently for subsequent discussion.

Proof Under the initial conditions (2.2), the existence and uniqueness of solutions of system (2.3) are to be proved step by step.

For 0 < t ≤ τ , it can be seen that I (t − τ ) ≥ 0 from the given initial condition. Then, we have that the right-hand side of system (2.3) is locally Lipschitz continuous. Therefore, by the existence, uniqueness and continuation theorems of differential equation, there is a unique solution of system (2.3) with the initial conditions (2.2) 

In what follows, we prove the nonnegativity of this solution in (0, τ ]. Taking notice of

> 0, and the initial condition S(ϑ) ≥ 0, we can obtain S(t) ≥ 0 for 0 < t ≤ τ . In the same way,

If not, we can find the smallest t 0 ∈ (0, τ ] which makes I (t 0 ) = 0 and I (t) < 0 for t ∈ (t 0 , t 0 + δ 1 ) where δ 1 > 0. In this case, together with the initial value I (ϑ) ≥ 0, we know I (t 0 −τ ) ≥ 0. Thus, there exists

Apparently, there is a contradiction in the interval (t 0 , t 0 + min{δ 1 , δ 2 }). Therefore, we draw a conclusion I (t) ≥ 0 for 0 < t ≤ τ . Utilizing the nonnegativity of I (t), the initial condition R(ϑ) ≥ 0 and

For τ < t ≤ 2τ , it can be known that I (t − τ ) ≥ 0 from the above discussion. Then, the existence and uniqueness of solutions can also be guaranteed in the interval (τ, 2τ ]. Utilizing the similar method of proof, we are able to obtain the nonnegativity of S(t), H (t),

This process can proceed if we adopt the same way in (2τ, 3τ ], (3τ, 4τ ] and so on. To sum up, it's proved that a unique solution (S(t), H (t), I (t), R(t)) of system (2.3) can continuously exist in the maximal interval [0, ∞).

+ is as follows

Proof Based on the proof of Lemma 1, it's obtained that S(t), H (t), I (t) and R(t) are nonnegative. As a result,

Observe that the sum of all equations of system (2.3) yields

(2.5)

Thus, we can obtain

Moreover, from the first equation of system (2.3), we gain

By the comparison principle, we have

In this manner, it's also acquired that

By equating the right sides of dS(t) dt , dH (t) dt and d I (t) dt to zero, we can easily verify that system (2.4) has a disease-free equilibrium point of the form

. This is regarded as an idealization where the disease completely disappear in the network. To find the threshold of disease propagation, we are to figure out the basic reproduction number which indicates the scale of the new infected resulted from an infected individual by contact. Based on the method of next-generation matrix [57] , define the new infection matrix F and the transition matrix V , namely

Hence, we can calculate the spectral radius of F V −1 and define the basic reproduction number of the infection as

Next, the condition for the existence of endemic equilibrium point E * of system (2.4) is determined in the following conclusion.

= 0 and suppose E * = (S * , H * , I * ) as the endemic equilibrium point of system (2.4). Hence, we obtain a quadratic equation about I * which satisfies the form q 0 I * 2 + q 1 I * + q 2 = 0 where

Define g(I * ) = q 0 I * 2 + q 1 I * + q 2 and note g(0) = (b 1 +b 2 +μ)(γ +μ)μ 2 a 2 (1− R 1 0 ) < 0 in the condition of R 1 0 > 1. Furthermore, we figure out

From the above discussion about g(I * ), we conclude that the equation g(I * ) = 0 has a unique positive root

, I * .

As a matter of fact, stability behavior is an important feature of the epidemic model as it reveals the stable characteristic of dynamical system in the long term. In this section, an earnest attempt is made to study the local and global stabilities of two equilibrium points of our delayed system.

Proof To investigate the local stability of the diseasefree equilibrium point, system (2.4) needs to be linearized at E 0 and the linear system has the following form

(2.7)

The Jacobian matrix J (E 0 ) of system (2.4) at E 0 can be written as

The equation (2.8) has two negative real roots equal to −μ. Besides, the remaining eigenvalue of matrix

In what follows, we discuss the characteristic roots of equation (2.8) for τ = 0 and τ > 0 separately.

(I) When τ = 0, all solutions of the characteristic equation (2.8) have negative real part if R 1 0 < 1. In this case, the disease-free equilibrium point E 0 of system (2.4) is locally asymptotically stable. On the contrary, there exists one eigenvalue with positive real part, which implies that E 0 is unstable in the case of R 1 0 > 1.

(II) For ∀τ > 0, assume that R 1 0 < 1 first. When instability at E 0 happens, one eigenvalue of the characteristic equation (2.8) is bound to cross the imaginary axis. Next, we wonder whether there is a pair of complex conjugate roots that passing through the imaginary axis with the increase of τ . Thus, suppose that there is a pair of purely imaginary roots. Substitute λ = iξ (ξ > 0) into the equation H (λ) = 0. By separating the real and imaginary parts, we obtain

(2.10)

Squaring and adding both sides of the two equations (2.10), we gain the following equation

In fact, there are no positive real roots ξ satisfying the above equation (2.11) when R 1 0 < 1.

Then, under the assumption of R 1 0 > 1, we note H (0) < 0 and figure out

Accordingly, find lim λ→+∞ H (λ) = +∞. As a result, the equation H (λ) = 0 has at least one positive real root if R 1 0 > 1, which suggests the instability at E 0 of system (2.4).

To summarize, combining the above two situations τ = 0 and τ > 0, we prove that the disease-free equilibrium point E 0 is locally asymptotically stable if

Proof According to Lemma 2, we find that all solutions of system (2.1) will remain or tend to the invariant region Ω 1 . In consequence, construct the following Lyapunov function L 1 (t) in the closed set Ω 1

In fact, we are to show that the derivative of L 1 (t) along the solutions of system (2.1) isn't positive for all t ≥ 0. Now, calculate the derivative of L 11 (t) with respect to t as follows:

Then, calculate the derivative of L 12 (t) with respect to t as follows:

Thus, the derivative of L 1 (t) along the solutions of system (2.1) is given by

When the basic reproduction number satisfies R 1 0 ≤ 1, we can obtain dL 1 

By calculation, the characteristic equation of J (E * ) is equivalent to

where

(I) When τ = 0, the characteristic equation (2.14) at E * can be degenerated into the form

Pay attention to the identical equation

(2.16) Therefore, we can find that coefficients c 1 , c 2 and c 3 of the equation (2.15) are all positive. Then, make some calculations in the following

According to the Hurwitz criterion, it can be proved that the endemic equilibrium point E * of system (2.4) is locally asymptotically stable in the case of τ = 0 when R 1 0 > 1. (II) When τ > 0, we make efforts to study the influence of time delay on the stability of system (2.4) at E * . Once system (2.4) generates instability for a specific delay τ , there exists a characteristic root of equation (2.14) which must cross the imaginary axis. Assume that λ = iξ (ξ > 0) is a solution of the equation (2.14), which meets the following form

(2.17)

By separating the real and imaginary parts of the equation (2.17), we derive n 5 ξ cos(ξ τ ) + (n 4 ξ 2 − n 6 ) sin(ξ τ ) = ξ 3 − n 2 ξ, (n 6 − n 4 ξ 2 ) cos(ξ τ ) + n 5 ξ sin(ξ τ ) = n 1 ξ 2 − n 3 .

After eliminating τ in (2.18), it follows that

Letting z = ξ 2 , we can rewrite (2.19) as an equation about z in the following z 3 + n 21 z 2 + n 22 z + n 23 = 0, (2.20) where n 21 = n 2 1 − 2n 2 − n 2 4 , n 22 = n 2 2 − 2n 1 n 3 − n 2 5 + 2n 4 n 6 and n 23 = n 2 3 − n 2 6 . Further, we are devoted to analyzing the property of coefficients n 21 , n 22 and n 23 . It's easy to find n 23 = (n 3 + n 6 )(n 3 − n 6 ) > 0 because n 3 + n 6 = β 1 k I * a + I * + μ β 2 k I * a + I * + μ (γ + μ)

Together with the identical equation (2.16), it's available that

And

Therefore, we prove that n 21 , n 22 and n 23 are all positive. Define h(z) = z 3 + n 21 z 2 + n 22 z + n 23 and note h(z) > 0 for ∀z > 0. The equation (2.20) has no positive roots z, which implies that there are no real roots ξ satisfying the equation (2.17). As a result, it can be concluded that the endemic equilibrium point E * of system (2.4) is locally asymptotically stable for ∀τ > 0 in the condition of R 1 0 > 1.

In this section, when the network is from homogeneous to heterogeneous, we propose and analyze the corresponding epidemic model which possesses the same modeling mechanism as that in the above section.

Compared with homogeneous networks, the network with a heterogeneous topology is closer to social network in the real world accurately, which motivates us to take the heterogeneity of contact between individuals into account. As a result, we make efforts to study disease transmission in heterogeneous networks in which nodes with the same degree are regarded as equivalents. Suppose that the network which we consider is uncorrelated. Denote the densities of the subhealthy, healthy, infected and recovered nodes with degree k (k = 1, 2, . . . , n) at time t as S k (t), H k (t), I k (t) and R k (t), respectively. The SHIR epidemic model in heterogeneous networks is under the framework of the following coupled system:

represents the density of newborn susceptible individuals coming from empty nodes. The definitions of all parameters in system (3.1) are same with that in systems (1.1) and (2.1). Most important of all, by introducing the psychological factor a, we significantly adopt a nonlinear incidence rate β 1 kS k (t) Θ(t) a+Θ(t) or β 2 k H k (t) Θ(t) a+Θ(t) to characterize the validity of protective measures taken by the subhealthy or healthy individuals. Furthermore, two susceptible groups need some time to be infected in the spreading process of disease. The consideration of delayed nonlinear incidence rate makes this epidemic model more reasonable and practical. The initial conditions (S k (ϑ), H k (ϑ), I k (ϑ), R k (ϑ)) k = (ϕ k1 (ϑ), ϕ k2 (ϑ), ϕ k3 (ϑ), ϕ k4 (ϑ)) k for system (3.1) satisfy

x k3 , x k4 ) k ∈ R 4n : x ki ≥ 0, i = 1, 2, 3, 4 .

Taking the sum of all 4n equations in system (3.1), we

Further, define T * k = lim t→∞ T k (t) = b 1 +b 2 b 1 +b 2 +μ . Then, system (3.1) can be transformed into the limit system as follows:

It can be noticed that the variable R k (t) doesn't affect values of S k (t), H k (t) and I k (t) in the above system. Moreover, the limit system (3.3) is further reduced to the following form:

(3.4)

In consequence, the subsystem (3.4) of the limit system (3.3) facilitates us to investigate dynamical behavior of system (3.1) sufficiently.

Proof Here, we intend to omit the specific proof which is similar to that of Lemma 1.

+ is as follows

Proof The detailed proof is omitted since it has a certain similarity to that of Lemma 2.

There is a unique endemic equilibrium point E * = (S * k , H * k , I * k ) k of system (3.4) when the basic reproduction number R 2 0 > 1.

Proof It's easily observed that system (3.4) always has a disease-free equilibrium point

which implies all of infected and recovered individuals finally degenerate into the susceptible state including the subhealthy and healthy. In addition, the endemic equilibrium point E * satisfies the following equations

. Further, the positive solution of the above equations is calculated as

,

,

.

Inserting the expression of I * k (t) into Θ * , we derive a self-consistent equation about Θ * , namely

The above equation (3.6) can be transformed into

By calculation, we have

Besides, there exists the following inequality

In addition, we figure out

Since system parameters are positive, it's verified that

To ensure that f (Θ * ) = 0 has a non-trivial solution Θ * ∈ (0, 1) , f (Θ * ) needs to satisfy the inequation f (0) > 0. Therefore, the basic reproduction number of the infection is expressed by

In conclusion, only if the condition R 2 0 > 1 is satisfied, system (3.4) admits a unique endemic equilibrium point E * = (S * k , H * k , I * k ) k in heterogeneous networks.

Theorem 6 For any τ ≥ 0, the disease-free equilibrium point E 0 of system (3.4) is locally asymptotically stable if R 2 0 < 1 and unstable if R 2 0 > 1.

Proof To study the local stability of system (3.4) at E 0 , we are aimed at deriving the linear system which is a good approximation to the nonlinear system (3.4). The linear dynamical system is expressed as

(3.9)

Correspondingly, the Jacobian matrix J (E 0 ) is a 3n × 3n matrix with the form of

and J 33 has the following form

Then, we can obtain the characteristic equation of the above Jacobian matrix J (E 0 )

(3.10)

As mentioned above, the characteristic equation (3.10) has 2n eigenvalues equal to −μ. In addition, n − 1 negative roots of the equation (3.10) are calculated as −(γ + μ). The remaining eigenvalue of matrix J (E 0 ) satisfies the following equation

(I) In the case of R 2 0 < 1, we wonder whether the eigenvalue λ can get to the imaginary axis in the plane or not. Suppose first that Reλ ≥ 0. From the equation (3.11), we take advantage of Euler's formula and derive the following result Reλ = (γ + μ) R 2 0 e −Reλτ cos(I mλτ ) − 1 < (γ + μ)(R 2 0 − 1) < 0, which contradicts the assumption apparently. Now, it's evidenced that all roots of the characteristic equation (3.10) have negative real part. For any τ ≥ 0, we demonstrate that the disease-free equilibrium point E 0 of system (3.4) is locally asymptotically stable if R 2 0 < 1.

(II) In the case of R 2 0 > 1, it's easy to make a calculation dG(λ) dλ

Taking notice of G(0) = −(γ + μ)(R 2 0 − 1) < 0 and lim λ→+∞ G(λ) = +∞, we can infer that there is at least one positive real root of the equation G(λ) = 0.

From the above discussion, it can be concluded that the disease-free equilibrium point E 0 is locally asymptotically stable for ∀τ ≥ 0 if R 2 0 < 1 and unstable if R 2 0 > 1. Theorem 7 For any τ ≥ 0, the disease-free equilibrium point E 0 of system (3.1) is globally asymptotically stable when R 2 0 ≤ 1. Proof In the invariant set Ω 2 , the Lyapunov function L 2 (t) is constructed as (3.12) where 

Similarly, we have

Hence, the derivative of L 2 (t) along the trajectories of system (3.1) is as follows

Obviously, it's available that dL 2 

is the singleton {E 0 }. When R 2 0 ≤ 1, the global asymptotic stability of disease-free equilibrium point E 0 of system (3.1) is proved according to the LaSalle Invariance Principle. That is to say, all solution trajectories initiating from different initial conditions arrive at the disease-free equilibrium point E 0 in the end, which implies that the infectious disease ultimately dies out regardless of the initial density of infected individuals in heterogeneous networks.

Remark 1 From (2.6) and (3.8), it's observed that the basic reproduction numbers R 1 0 and R 2 0 depend on the birth (death) rates, transmission coefficients, recovery rate, psychological factor and network structure. Comparing the conditions for the global stability of diseasefree equilibrium points E 0 and E 0 , we find that R 2 0 ≤ 1 is harder to be satisfied than R 1 0 ≤ 1 under the same parameters.

In this section, numerical simulations are carried out to illustrate analytical results obtained above. Besides, we also attempt to complement dynamics of the SHIR epidemic model with time delay. More precisely, system (2.4) is used as the object for simulations in homogeneous networks. For system (3.4) , the scale-free network with p(k) ∝ 2m 2 k −3 is considered as the underlying network. Here, the minimum and maximum degrees of nodes are, respectively, assumed to be m = 1 and n = 200. Moreover, the average degree of all nodes is computed as k = 1.9802 in the network.

As a crucial threshold, the basic reproduction number can determine the existence and stability of equilibrium points of dynamical systems. Hence, we are to observe the influence of model parameters on the basic reproduction number. Firstly, we focus on presenting the influence of average degree k and psychological factor a on R 1 0 in homogeneous networks together. Model parameters are chosen as b 1 = 0.06, b 2 = 0.09, β 1 = 0.35, β 2 = 0.28, γ = 0.8, μ = 0.5, a ∈ [0.1, 1] and k ∈ [0.1, 60]. As Fig. 1 displays, different color codings characterize different values of basic reproduction number R 1 0 vividly. It can be observed that R 1 0 increases as k increases or a decreases. This result suggests that the epidemic situation becomes more severe when an individual has the greater average number of neighbors in the social network. Nevertheless, good psychological quality of the public facing the epidemic can be beneficial to reduce the transmission of infectious disease. In addition, when the transmission rate β 1 or β 2 gets larger, the basic reproduction number R 2 0 becomes larger, which implies reducing the transmission rate from the subhealthy or healthy to infected individuals can weaken the spread of disease. With the increase of psychological factor, recovery or death rate, the basic reproduction number R 2 0 decreases gradually. Besides, the heterogeneity of underlying network expressed by k 2 /k contributes to R 2 0 , which reveals the decline of heterogeneity is conductive to eradicating the disease propagation. The influence of network heterogeneity on R 2 0 can also be observed by replacing k with k 2 /k in the abscissa of Fig. 1. In Fig. 2a , set the following parameters β 1 = 0.5, β 2 = 0.05, a = 0.86, γ = 0.77, μ = 0.55, b 1 ∈ [0, 1] and b 2 ∈ [0, 1]. Under the given parameters, the larger b 2 is, the smaller R 2 0 is. Similarly, the larger β 1 or b 1 can result in the growth of R 2 0 , as shown in Fig. 2b , which implies a higher birth rate of subhealthy individuals makes the disease easier to spread. Figure 3 depicts the changing trend of mean densities of subhealthy and infected nodes with the elapse of time. Likewise, take parameters b 1 = 0.5, b 2 = 0.53, β 1 = 0.6, β 2 = 0.3, a = 0.9, γ = 0.32, μ = 0.42, τ = 2 and same initial conditions. The basic reproduction number is given by R 1 0 = 0.9412 < 1. Figure 3 shows that the mean density of subhealthy nodes stabilizes to a constant state and that of infected nodes reduces to zero gradually when R 1 0 < 1, which implies the disease wipes out at last. From Theorem 2, the disease-free equilibrium point E 0 is locally asymptotically stable, which is in line with the observed result in Fig. 3 to some extent.

To verify the global asymptotic stability of the disease-free equilibrium point E 0 , we would like to observe the influence of initial conditions on the spread of disease by choosing 13 sets of different initial densities of subhealthy, healthy and infected individuals. Besides, model parameters are taken as b 1 = 0.4, b 2 = 0.5, β 1 = 0.82, β 2 = 0.5, a = 0.6, γ = 0.92, μ = 0.55, τ = 3.5. Under this condition, the basic reproduction number is R 1 0 = 0.8950 < 1 and the disease-free equilibrium point of system (2.4) is E 0 = (0.276, 0.345, 0) . What's more, it can be seen in Fig. 4 that all trajectories of system (2.4) in homogeneous networks ultimately converge toward E 0 , which is consistent with Theorem 3. nodes are to keep on the respective positive levels when R 1 0 > 1, which illustrates the disease is persistent in homogeneous networks. In addition, we discover the mean density of S(t) becomes smaller while that of I (t) gets bigger with R 1 0 increasing. In other words, the larger R 1 0 gives rise to the higher level of disease propagation, which suggests we can control the spread of disease effectively by decreasing the basic reproduction number. 4) . The basic reproduction number is calculated as R 2 0 = 3.5103 > 1. Figure 11 portrays evolutions of healthy and infected nodes with different k = 20, 60, 100, 140 and 180, respectively. With the rise of k, the density of healthy individuals decreases whereas the density of infected individuals increases. In reality, this phenomenon shown in Fig. 11 is actually understandable since nodes with more neighbors have easier accesses to be infected with the disease by contacting with the infected. In heterogeneous networks, regarded as hot spots, nodes connected with quantities of edges have major roles in propagating the disease. This result provides us with the idea to focus on this 4) . Then, the basic reproduction number is R 1 0 = 6.6007 > 1 correspondingly. By choosing τ = 0 and 10, Fig. 12 presents the difference between the evolutions of individuals without and with time delay. It can be observed that time delay τ has a remarkable influence on the convergence rates of S(t), H (t) and I (t). To be specific, once hysteresis exists in the infection process, more time is needed for system (2.4) to approach the endemic steady state. As also shown in Fig. 12 , the peak of the density of infected individuals I (t) with delay is lower than that without delay, which implies that infection delay weakens the maximum prevalence of disease. Besides, set b 1 = 0.24, b 2 = 0.36, β 1 = 0.31, β 2 = 0.25, a = 0.9, γ = 0.86, μ = 0.5 and S k (0) = 0.15, H k (0) = 0.1, I k (0) = 0.2 in system (3.4) . Under this condition, the basic reproduction number is R 2 0 = 0.9602 < 1. For different τ = 0.5, 1.5, 2.5, 3.5 and 4.5, densities of healthy and infected individuals, respectively, converge to the steady state H 0 k = 0.327 and I 0 k = 0 at last in Fig.  13 . Nevertheless, system needs more time to gather to the disease-free equilibrium point E 0 with the increase of τ . As a result, diminishing the time delay in the process of disease propagation can make the epidemic enter a stable situation at a faster speed. Through the simulation, it's shown that systems (2.4) and (3.4) produce only stable steady-state dynamics, which implies that the time delay incorporated in nonlinear incidence rate is harmless in our model essentially. S(t) with bilinear incidence rate H(t) with bilinear incidence rate I(t) with bilinear incidence rate S(t) with nonlinear incidence rate H(t) with nonlinear incidence rate I(t) with nonlinear incidence rate With the same parameters, there exists an obvious difference in the value of basic reproduction number, which can be attributed to the complexity of network structure. In heterogeneous networks, the total density of subhealthy nodes at time t is defined as S(t) = n k=1 p(k)S k (t). Similarly, H (t) = n k=1 p(k)H k (t) and I (t) = n k=1 p(k)I k (t) denote the total densities of healthy and infected nodes at time t. As displayed in Fig. 16, S(t) , H (t) and I (t) converge more slowly in homogeneous networks than that in heterogeneous networks due to R 1 0 < R 2 0 . Furthermore, the density of infected nodes tends to zero in homogeneous networks while that approaches the endemic steady state in heterogeneous networks. Namely, even if system parameters and initial conditions are same, the disease disappears eventually in 16 Evolutions of subhealthy, healthy and infected nodes in both homogeneous and heterogeneous networks homogeneous networks, whereas it still prevails in heterogeneous networks. This reveals that the heterogeneity of network is capable to exacerbate the propagation of disease. If the heterogeneous network is divided into several homogeneous networks, the prevalence of disease will be reduced to some extent.

Based on an SIR model with two susceptible groups in Ref. [30] , in this paper we come up with a modified SHIR epidemic model with time delay and nonlinear incidence rate in complex networks. In virtue of the compartment approach, two susceptible groups are interpreted as subhealthy and healthy individuals in view of the difference in fitness levels which indirectly affects the susceptibility of individuals to infectious disease. By introducing the psychological factor a, we adopt the nonlinear incidence rate to reflect the behavioral changes of susceptible individuals due to the psychological effect during the progression of disease diffusion. Besides, time delay τ is taken into account in the transformation process from the susceptible to the infected, which poses a challenge to the stability analysis of our network-based epidemic model. Specifically, it's logical for us to think about two scenarios incorporating τ = 0 and τ > 0 in studying the local asymptotic stability of equilibrium points. Meanwhile, the construction of Lyapunov function is also closely related to the time delay when analyzing the global asymptotic stability of equilibrium points. In addition, the topological structure of network is considered in homogeneous and heterogeneous networks, respectively. Some main results of the delayed SHIR epidemic model in complex networks are presented as follows.

(1) In homogeneous networks, we obtain equilibrium points and a positive invariant set of system based on the mean field equations. The basic reproduction number R 1 0 of the infection is defined according to the method of next-generation matrix. In the condition of R 1 0 < 1, we prove that the disease-free equilibrium point E 0 is both locally and globally asymptotically stable. Furthermore, the proof of the local asymptotic stability of E * is given in detail when R 1 0 > 1. (2) In heterogeneous networks, for the subsystem of the limit system, we investigate the disease-free equilibrium point E 0 = The basic reproduction number R 2 0 is derived on the basis of the existence of a positive equilibrium point. By means of linearization, we demonstrate the local asymptotic stability of E 0 if R 2 0 < 1 and instability of E 0 if R 2 0 > 1. Additionally, by constructing suitable Lyapunov function, the global asymptotic stability of E 0 is also studied when R 2 0 ≤ 1. (3) To examine and complement analytical results, quantities of numerical simulations are executed orderly. We conduct the sensitivity analysis of model parameters on the basic reproduction number. In particular, the heterogeneity of network structure makes signifi-cant difference in disease propagation. Then, stability of equilibrium points in complex networks is verified intuitively. Further, we devise several numerical experiments to explore the potential dynamical behavior at the endemic equilibrium points in both homogeneous and heterogeneous networks. Moreover, the influence of key parameters on disease diffusion mainly reflects in the speed of propagation and steady-state densities of individuals.

In fact, the study of our two-susceptibility epidemic model with time delay and nonlinear incidence rate lays the foundation for studying the more complicated dynamical system with multiple groups later. Driven by the practical significance of epidemic models, we hope to access actual data to carry out statistical research in the future.

Analysis of a reactiondiffusion cholera epidemic model in a spatially heterogeneous environment

Modeling and analyzing dynamics of malaria transmission with host immunity

A fractional order epidemic model and simulation for avian influenza dynamics

The reproductive number of COVID-19 is higher compared to SARS coronavirus

Stability analysis of an SIS epidemic model with feedback mechanism on networks

Dynamics of epidemic spreading model with drug-resistant variation on scale-free networks

Dynamic analysis of a delayed model for vector-borne diseases on bipartite networks

Epidemic spreading and global stability of an SIS model with an infective vector on complex networks

Stability analysis and control strategies for worm attack in mobile networks via a VEIQS propagation model

The dynamics analysis of a rumor propagation model in online social networks

Partial differential equation modeling of rumor propagation in complex networks with higher order of organization

A novel dynamic model for web malware spreading over scale-free networks

Web malware spread modelling and optimal control strategies

Modeling of knowledge transmission by considering the level of forgetfulness in complex networks

Knowledge transmission model with consideration of self-learning mechanism in complex networks

Dynamical behavior of a stochastic multigroup SIR epidemic model

An intuitive formulation for the reproductive number for the spread of diseases in heterogeneous populations

Differential susceptibility epidemic models

Global stability of an SIR model with differential infectivity on complex networks

The spreading dynamics of sexually transmitted diseases with birth and death on heterogeneous networks

Further dynamic analysis for a network sexually transmitted disease model with birth and death

Transmission dynamics of a multi-group brucellosis model with mixed cross infection in public farm

Modeling direct and indirect disease transmission using multi-group model

The dynamics of spreading and immune strategies of sexually transmitted diseases on scale-free network

Dynamics of a timeperiodic two-strain SIS epidemic model with diffusion and latent period

Analytical and numerical approaches to coexistence of strains in a two-strain SIS model with diffusion

Global results for an HIV/AIDS model with multiple susceptible classes and nonlinear incidence

Global stability of a multigroup SEI model

A deterministic model for gonorrhea in a nonhomogeneous population

Global stability of an SIR model with two susceptible groups on complex networks

Dynamics of a delay differential equation model of hepatitis B virus infection

Mathematical analysis of the global dynamics of an eco-epidemiological model with time delay

Global stability analysis of an SIR epidemic model with demographics and time delay on networks

Global dynamics analysis and control of a rumor spreading model in online social networks

Super-critical and critical traveling waves in a two-component lattice dynamical model with discrete delay

A delayed avian influenza model with avian slaughter: stability analysis and optimal control

Dynamics for an SEIRS epidemic model with time delay on a scale-free network

Delay differential equations modeling of rumor propagation in both homogeneous and heterogeneous networks with a forced silence function

Dynamical behaviours and control measures of rumour-spreading model with consideration of network topology

Nonlinear dynamics of a timedelayed epidemic model with two explicit aware classes, saturated incidences, and treatment

Endemic bubbles generated by delayed behavioral response: global stability and bifurcation switches in an SIS model

A generalization of the Kermack-McKendrick deterministic epidemic model

Global analysis of an epidemic model with nonmonotone incidence rate

Dynamics of a network-based SIS epidemic model with nonmonotone incidence rate

Dynamics of a time-delayed SIR epidemic model with logistic growth and saturated treatment

Bifurcation analysis of an SIV epidemic model with the saturated incidence rate

A deterministic time-delayed SIR epidemic model: mathematical modeling and analysis

Nonlinear dynamical analysis and control strategies of a network-based SIS epidemic model with time delay

Bifurcation analysis of a networkbased SIR epidemic model with saturated treatment function

Stability and extinction of SEIR epidemic models with generalized nonlinear incidence

Rumor spreading model considering hesitating mechanism in complex social networks

Rumor spreading model with noise interference in complex social networks

Dynamical analysis of a rumor spreading model with self-discrimination and time delay in complex networks

Review mechanism promotes knowledge transmission in complex networks

Stability analysis of I2S2R rumor spreading model in complex networks

Stability behavior of a nonlinear mathematical epidemic transmission model with time delay. Nonlinear Dyn

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

Data availability This manuscript has no associated data.

The authors declare that there is no conflict of interests regarding the publication of this article.

The authors state that this article complies with ethical standards. This article does not contain any studies with human participants or animals.