Complex dynamics of an epidemic model with saturated media coverage and recovery.

During the outbreak of emerging infectious diseases, media coverage and medical resource play important roles in affecting the disease transmission. To investigate the effects of the saturation of media coverage and limited medical resources, we proposed a mathematical model with extra compartment of media coverage and two nonlinear functions. We theoretically and numerically investigate the dynamics of the proposed model. Given great difficulties caused by high nonlinearity in theoretical analysis, we separately considered subsystems with only nonlinear recovery or with only saturated media impact. For the model with only nonlinear recovery, we theoretically showed that backward bifurcation can occur and multiple equilibria may coexist under certain conditions in this case. Numerical simulations reveal the rich dynamic behaviors, including forward-backward bifurcation, Hopf bifurcation, saddle-node bifurcation, homoclinic bifurcation and unstable limit cycle. So the limitation of medical resources induces rich dynamics and causes much difficulties in eliminating the infectious diseases. We then investigated the dynamics of the system with only saturated media impact and concluded that saturated media impact hardly induces the complicated dynamics. Further, we parameterized the proposed model on the basis of the COVID-19 case data in mainland China and data related to news items, and estimated the basic reproduction number to be 2.86. Sensitivity analyses were carried out to quantify the relative importance of parameters in determining the cumulative number of infected individuals at the end of the first month of the outbreak. Combining with numerical analyses, we suggested that providing adequate medical resources and improving media response to infection or individuals' response to mass media may reduce the cumulative number of the infected individuals, which mitigates the transmission dynamics during the early stage of the COVID-19 pandemic.

Introduction

The emerging or re-emergent infectious diseases, including SARS(2003), HINI(2009) and COVID-19, have been threatening the public health and caused worrying concern amongst the public and health authorities. In the initial stage of the epidemic, massive news coverage and fast information flow significantly generate profound psychological/behavioural impacts on the public, which may influence disease spread and then the implementation of public interventions [1]. There are a number of mathematical models considering the media impact on disease. Some included a factor in the incidence rate to lower the infection [2-6]. Some proposed the on-off media functions [5-7], where the media impact is triggered by the disease infection, they show that media coverage can be considered as an effective way to mitigate the disease spreading during the initial stage of an outbreak. Others modelled the level of media coverage as an independent compartment [4, 8–12]. Song et al. [9] further extended the models by considering the delayed media function and examined local/global bifurcation caused by time delays. Wang et al. [5] proposed an epidemic model with the on-off media function triggered by number of infected individuals, and testified that a previously chosen level of the desired number of infected individuals can be reached if the threshold policy and other parameters are chosen properly. Recently, Tiwari et al. [13] proposed a model which portrays the interplay between dissemination of awareness at local and global levels, and prevalence of disease, where global awareness stem from social media advertisements about disease. Many complex dynamic phenomena can occur in such model. It testified that global awareness is more effective in curtailing the disease than local awareness.

The limitation of medical resources is the most serious threat during the emerging infectious diseases, and occurs in almost all countries. For example, with sharp the increase of COVID-19 cases in mainland China in early 2020 the shortage of testing kits and hospital beds was the most serious problem. The influence of medical resources on the spread of infectious diseases has been paid an increasing attention. Cui et al. [14, 15] formulated an epidemic model with saturation recovery from infective individuals and showed that saturation recovery leaded to bistability and periodicity. Recently, to analyze the effect of medical conditions, Li and Zhang [16] proposed a SIR model in epidemic diseases using nonmonotone incidence and saturated recovery rates and investigated complex dynamics such as the saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation of codimension 2. It showed that a sufficient number of the beds is critical to control the epidemic. Asamoah et al. [17] proposed a mathematical model with nonlinear recovery rate for bacterial meningitis to provide a potential framework for the control of the disease spread in limited-resource settings, and found the forward and backward bifurcation. This indicated that nonlinear recovery (and/or nonmonotone incidence) may induce complicated dynamics.

A common assumption in above studies about media coverage is that the growth of media items linearly depends on the number of infected individuals. However, this is not how the thing looks like. It is known that, facing the emerging or re-emerging infectious diseases, news items surge in a certain period but have saturation due to the fatigue of reporting, and consequently the saturated growth of news items should be considered especially during whole epidemic outbreak. For example, with sharp the increase of COVID-19 cases in mainland China [18, 19], information dissemination significantly affected on disease spread via profound psychological/behavioural changes in early 2020, while it has been exhibiting the reporting fatigue recently. Further, with the media impacts, whether the nonlinear recovery induces more complex dynamics, and the detailed contributions of saturated media impact and/or nonlinear recovery to rich dynamics or disease control remain unclear, and hence fall within the scope of this study.

Our main purpose of this study is to propose a mathematical model with double nonlinear functions to investigate the impacts of saturated media coverage and limited medical resources on disease transmission and examine the rich dynamics induced by the nonlinearity. We initially consider model with only nonlinear recovery in Sect. 3. Theoretical analysis shows that backward bifurcation occurs under certain parameters, and multiple equilibria can coexist in such system and numerical simulations reveal some rich dynamic behaviors, including forward-backward bifurcation, Hopf bifurcation, saddle-node bifurcation, homoclinic bifurcation and unstable limit cycle. Then we analyze model with only saturated media coverage in Sect. 4 to clear the contribution of saturated media on dynamics, in which the threshold dynamics are investigated and the global stability of the endemic state under certain conditions is proved. In Sect. 5, we parameterize the proposed model on the basis of the COVID-19 case data in mainland China and the data related to news items, and estimate the basic reproduction number. Sensitivity analysis is carried out to quantify the relative importance of parameters in determining the cumulative number of newly observed infectious individuals at the end of the first month of outbreak. Combining with numerical simulations, we testify that the effectiveness of providing adequate medical resources and improving media response to infection or individuals’ response to mass media. Finally, a brief discussion is given.

Model description

We divide the population into susceptible individuals (S), exposed but not yet infectious individuals (E) and infected individuals (I). We further assume that the susceptible individuals are infected by infectious individuals with a rate of , and become exposed. The exposed individuals become infectious with a rate of , and then recover with a saturated function and back into susceptibles class, where is the recovery rate of the infective individuals and is the parameter that measures the effect of medical resource limitation [14, 15].

In order to evaluate the effect of media coverage and fast information flow on disease transmission, we include the average number of daily news items related to the outbreak as an independent variable, denoted as M(t), the changing rate of which depends on the number of newly observed infectious individuals () with a saturated function , where represents the reporting rate and is the parameter that measures the saturation effect of media reports. We further assume that media reports wane at the rate of . To respond to the effect of media coverage on reducing the effective contact of susceptible with infectious individuals, we further include an exponential decreasing factor [9, 20], which is represented by , to describe the reduction factor in the incidence. Thus the model equations arewhere is the rate of flowing into the population and d is the natural death rate. All parameters are non-negative constants.

We can testify that the feasible region of model (2.1) isThe disease-free equilibrium is always feasible. By calculating the spectral radius of the next generation for the model (2.1), we get the basic reproduction number of this model iswhere is the expected number of secondary cases produced, in a completely susceptible population, by a typical infected individual during his/her infectious period [21]. The local stability of disease-free equilibrium also can be proved.

Theorem 1

The disease-free equilibrium is locally asymptotically stable if and unstable if . The bifurcation at is backward when ; the bifurcation at is forward when , that is to say when , is locally asymptotically stable if and unstable if , where

Proof

The Jacobian matrix concerned with the linearization of system (2.1) at isSo, the corresponding characteristic equation isThis implies that the disease-free equilibrium is locally asymptotically stable if and unstable if .

When , zero is a simple eigenvalue of and all other eigenvalues of have negative real parts. Hence is a nonhyperbolic equilibrium, and the linearization cannot determine its stability. Here, we use center manifold [22] to analyze its stability. Let , where is the transpose of the matrix, so model (2.1) can be denoted as . We denoted the left and right eigenvalues of as v and w respectively, where , . Algebraic calculations show thatwhere . One can further choose m and n such that and .

Let and , then if and only if . According to Theorem 4.1 in [23], the local dynamic behavior near is determined by the following constants a and b:In terms of a and b, we just need to find the nonzero partial derivative of and at when and . Algebraic calculations show thatAll other partial derivative of and are zero. So, we can get thatFurther, when , must hold, then we can get if and only ifSo we get the bifurcation at is backward when ; the bifurcation at is forward when , that is to say when , is locally asymptotically stable if and unstable if . In particular, when , holds true, so the bifurcation at is forward and backward bifurcation can not occur in such scenario. That completes the proof.

It is worth noticing that the system (2.1) may occur backward bifurcation at , which is mainly associated with parameter , instead of parameter . This means that nonlinear recovery function significantly contributes the occurrence of backward bifurcation. We can also easily prove the persistence of the system (2.1) for by using the standard methods [24], shown in the following theorem (see Appendix A.1 for details).

Theorem 2

If , there exists a positive constant such that for the system (2.1)

Note that system (2.1) has high nonlinearity, which is complicated to analyze the dynamical behaviors of this system and the occurring of backward bifurcation for system (2.1) at is mainly associated with parameter , instead of parameter , we initially consider a special case: model with only nonlinear recovery (i.e., ), and conclude the dynamics of the proposed system theoretically and numerically.

Dynamical behaviors of the system (2.1) with

When only considering the saturated recovery (i.e., ), the corresponding model equations areThe disease-free equilibrium of system (3.5) is the same as that for system (2.1), and in the following we focus on the stability of the DFE, the feasibility of the endemic states and possible bifurcations around them.

Existence of the endemic equilibrium

As for system (3.5), the possible endemic equilibrium satisfiesDenote , andThen the existence of the solution of the last equation in (3.6) can be translated into the existence of the intersection of the functions and . Algebraic calculation yieldsfor , whereFurther, we can testify thatif and only if andwhere .

When , that is to saywe have for , so must holds for in this case, which implies that function decreases monotonically on for . Considering the zero point theorem and combining the monotonicity of for , we conclude that function and have a unique intersection on the interval for and ; function and have no any intersection on the interval for and .

When , that is to say . Denote , we get and for , then must hold true for , so we just need to consider if functions and have an intersection on . For , must hold, which implies that the system (3.5) at least has one endemic equilibrium for on the basis of the zero point theorem. For , must hold. Denote . When , functions and have no intersection on . For , there must be a , such that , so and . According to the zero point theorem, f(I) has at least two zero point, one is in the interval of and the other in , i.e., the system has at least two positive equilibriums.

According to the above analysis, we obtain the following conclusion.

Theorem 3

If , the system (3.5) has a unique endemic equilibrium for , denoted by , which satisfies (3.6) and ; but it has no endemic equilibrium for . If , the system (3.5) exists at least one endemic equilibrium for ; the system (3.5) has no endemic equilibrium for and , but it at least exists two endemic equilibrium for and .

Theorem 4

When and , the unique endemic equilibrium of system (3.5) is locally asymptotically stable.

Proof

The Jacobian matrix concerned the linearization of system (3.5) at isLet , we can getandThen the characteristic equation yieldswithwhere .

Further, we havewhere .

We can get the system (3.5) has a unique endemic equilibrium for in the above cases, on the basis of Theorem 3. And must be true.

When , the inequalityholds true sine the right side of inequality (3.13) is less than 0 if and the left side of inequality (3.13) must be greater than 0. And the inequality (3.13) is equivalent to the inequalitywhere . The inequality (3.14) is further equivalent to the following inequalitywhich is naturally satisfied. Then in this case, () and must hold true from the expressions of () and .

So if , () and hold true. According to the Hurwitz criterion, the roots of the polynomial all have negative real parts, which means that the roots of the characteristic equation (3.10) all have negative real parts. So this unique endemic equilibrium is locally asymptotically stable for . That completes the proof.

Theorem 5

The disease-free equilibrium of system (3.5) is globally asymptotically stable for and with .

Let’s consider a Lyapunov functionFor , we can obtainSo, if , we have , and holds if and only if . Applying LaSalle’s invariance principle, we get that . Then applying comparison principle, we can testify and . So the disease-free equilibrium is globally asymptotically stable if and . That completes the proof.

Note that system (3.5) has a unique positive equilibrium, which is locally stable, when and , we can further develop the global result of this system under certain parameter restrictions by using compound matrices and geometric ideas [25-28].

Theorem 6

When and , the unique endemic equilibrium of the system (3.5) is globally asymptotically stable if , and .

The Jacobian matrix of system (3.5) readsThe third additive compound matrix of J isand the corresponding linear compound system readswhereThen we just need to testify the global stability of this linear compound system according to [25, 26]. To do this, we choose the Lyapunov functional aswhere , will be determined later. It is easy to show the system (3.5) admits a unique nonnegative bounded solution for all and any non-negative initial value. Combing the Theorem 2, there exist two positive constants , such thatRecalling thatThen, discussing the derivative of V along the trajectory of the linear compound system (3.16) for different conditions and denoting the right-hand (total) derivative with respect to t.

Case (1): and , which implies . Algebraic calculation yieldsCase (2): and , which implies . Algebraic calculation yieldsSince , must hold, which means .

Case (3): and , which implies . Algebraic calculation yieldsCase (4): and , which implies . Algebraic calculation yieldsCase (5): and , which implies . Algebraic calculation yieldsSince , must hold, which means

Case (6): and , which implies }. Algebraic calculation yieldsSince , and , we can choose two positive constants and such that and respectively, which means that , , , and . Therefore, there exists a positive constant such that , in this case, must hold. By Corollary 3.2 in [26], of the system (3.5) is globally asymptotically stable. That completes the proof.

To further analyze the system (3.5) in detail, we investigate the number of positive equilibria and their stability for on the basis of Theorem 3. Here we have only analyzed a few specific cases, and others are similar but complicated.

Table 5

The concavity of

Cases	Conditions	Conclusion
Case A	(A.3) holds Q>0 and I~<Λd	f2^′′(I)>0 for 0<I<I~ and I~<I<Λd f2^′′(I)=0 for I=I~
Case B	(A.3) holds Q>0 and I~>Λd	f2^′′(I)>0 for 0<I<Λd
Case C	(A.5) holds Q>0 and I_2~<Λd	f2^′′(I)>0 for 0<I<I_1~ and I_2~<I<Λd f2^′′(I)<0 for I_1~<I<I_2~
Case D	(A.5) holds Q>0 and I_1~<Λd<I_2~	f2^′′(I)>0 for 0<I<I_1~ f2^′′(I)<0 for I_1~<I<Λd
Case E	(A.5) holds Q>0 and Λd<I_1~	f2^′′(I)>0 for 0<I<Λd
Case F	(A.6) holds and Q>0	f2^′′(I)>0 for 0<I<Λd
Case G	Q<0 and I~<Λd	f2^′′(I)<0 for 0<I<I~ f2^′′(I)>0 for I~<I<Λd
Case H	Q<0 and I~>Λd	f2^′′(I)<0 for 0<I<Λd

Let’s initially consider the case when and . Algebraic calculation yields the quadratic equation has a positive real root, denoted by , and a negative real root,where and satisfies . Then we have for and for . Further, we havewe then obtain that the function is a strictly convex function that increases on the interval and decreases on the interval . Similarly we can analyze the concavity of the function , which is given in Table 5 by considering eight cases (case A-H), see Appendix A.2 for details.

Combining the monotonicity and concavity of functions and , we discuss the number of local equilibria in the following three scenarios according to the number of intersection of functions and in the interval of (Table 1).

Table 1

The number of endemic equilibria

scenarios	Conditions	Concavity	Number of equilibria for R_0>1	Number of equilibria for R_0<1
scenario 1	Case B :(A.3) holds Q>0 and I~>Λd	f2^′′(I)>0 for 0<I<Λd	1	0 or 1 or 2
	Case E : (A.5) holds Q>0 and Λd<I_1~
	Case F : (A.6) holds and Q>0
	Case A : (A.3) holds Q>0 and I~<Λd	f2^′′(I)>0 for 0<I<I~ and I~<I<Λd f2^′′(I)=0 for I=I~
scenario 2	Case C : (A.5) holds Q>0 and I_2~<Λd	f2^′′(I)>0 for 0<I<I_1~ and I_2~<I<Λd f2^′′(I)<0 for I_1~<I<I_2~	1 or 2 or 3 (See Table 6 for details))	0 or 1 or 2 or 3 or 4
	Case D : (A.5) holds Q>0 and I_1~<Λd<I_2~	f2^′′(I)>0 for 0<I<I_1~ f2^′′(I)<0 for I_1~<I<Λd	1 or 2 or 3
scenario 3	Case G : Q<0 and I~<Λd	f2^′′(I)<0 for 0<I<I~ f2^′′(I)>0 for I~<I<Λd	1 or 2 or 3	0 or 1 or 2
	Case H : Q<0 and I~>Λd	f2^′′(I)<0 for 0<I<Λd

Case A, B, ... and H are given for discussing concavity of function , referring to Table 5

Scenario 1): If cases B, E or F hold, that is to say, is a strictly concave function of monotonically increasing. And the function is a strictly convex function that increases on the interval and decreases on the interval . Combining (3.7) and (3.8), we can conclude the system (3.5) has a unique endemic equilibrium for . When case A holds, is a concave function of monotonically increasing and we can do the similar analysis.

Scenario 2): If case C holds, system (3.5) may have 1 to 3 endemic equilibria we use a table (i.e., Table 6) to illustrate for simplicity. For the sake of easy to representation, if there is one, two or three endemic equilibria, we denote them as ; and (); , and () respectively. As for the case D, we can do the same analysis.

Table 6

The number of endemic equilibria of the case 2

Conditions	Amount	Range
I_1≤I_1~<I_2~	One	I^∗>I_1 for f_1(I_1)>f_2(I_1) I^∗<I_1 for f_1(I_1)<f_2(I_1) I^∗=I_1 for f_1(I_1)=f_2(I_1)
I_1~<I_1<I_2~, f_2(I_1~)>f_1(I_1~), f_1(I_1)>f_2(I_1)	Three	I1∗<I_1~<I2∗<I_1<I3∗<I_2~
I_1~<I_1<I_2~, f_2(I_1~)>f_1(I_1~), f_1(I_1)<f_2(I_1)	One	I^∗<I_1~
I_1~<I_1<I_2~, f_2(I_1~)>f_1(I_1~), f_1(I_1)=f_2(I_1)	Two	I1∗<I_1~<I2∗=I_1
I_1~<I_1<I_2~, f_2(I_1~)=f_1(I_1~), f_1(I_1)>f_2(I_1)	Two	I1∗=I_1~<I_1<I2∗
I_1~<I_1<I_2~, f_2(I_1~)=f_1(I_1~), f_1(I_1)<f_2(I_1)	One	I^∗=I_1~
I_1~<I_1<I_2~, f_2(I_1~)=f_1(I_1~), f_1(I_1)=f_2(I_1)	Two	I1∗=I_1~<I2∗=I_1
I_1~<I_1<I_2~, f_2(I_1~)<f_1(I_1~)	One	I^∗>I_1
I_1~<I_2~<I_1, f_2(I_1~)>f_1(I_1~), f_1(I_2~)>f_2(I_2~)	Three	I1∗<I_1~<I2∗<I_2~<I3∗
I_1~<I_2~<I_1, f_2(I_1~)>f_1(I_1~), f_1(I_2~)<f_2(I_2~)	One	I^∗<I_1~
I_1~<I_2~<I_1, f_2(I_1~)>f_1(I_1~), f_1(I_2~)=f_2(I_2~)	Two	I1∗<I_1~<I2∗=I_2~<I_1
I_1~<I_2~<I_1, f_2(I_1~)=f_1(I_1~), f_1(I_2~)>f_2(I_2~)	Two	I1∗=I_1~<I2∗
I_1~<I_2~<I_1, f_2(I_1~)=f_1(I_1~), f_1(I_2~)<f_2(I_2~)	One	I^∗=I_1~
I_1~<I_2~<I_1, f_2(I_1~)=f_1(I_1~), f_1(I_2~)=f_2(I_2~)	Two	I1∗=I_1~<I2∗=I_2~
I_1~<I_2~<I_1, f_2(I_1~)<f_1(I_1~)	One	I^∗>I_1

Scenario 3): If case G or H holds, we can get the system (3.5) has one or two or three endemic equilibria for similarly.

As for the case when , , we can use the same method to prove that system (3.5) may have 0 to 4 endemic equilibria (see details for Table 1), which means that the system (3.5) may occur backward bifurcation.

Due to computational difficulties, we will implement numerical simulations to analyze dynamical behaviors of (3.5) and possible bifurcations under various cases with the help of the MATCONT package [29]. We shall illustrate forward bifurcation, backward bifurcation, forward-backward bifurcations, Hopf bifurcation and saddle-node bifurcations in the following.

Fig. 4

a: Bifurcation diagram in () plane with , where , , , , , , , , . The blue curves represent stable equilibria, the red curves represent unstable equilibria and a black dot marked with “SN” ( and ) or “H” ( and ) indicates that the system (3.5) goes through saddle-node bifurcation or Hopf bifurcation. c and d: Local bifurcation diagram, where “Hom” ( and ) indicates that the system (3.5) goes through homoclinic bifurcation. d and e: The variation in period of cycles versus

A forward bifurcation is illustrated in Fig.1, and Fig.1(a) gives variation in equilibrium level of infections () with the basic reproduction number (associated with parameter ) and corresponding bifurcations, say, () and when (), the system (3.5) undergoes Hopf bifurcation. Figure 1(b-e) illustrates the time series of the system (3.5) with different parameters which stabilize to either the disease-free equilibrium or endemic state or bifurcated periodic solution. The first Lyapunov coefficient at two Hopf points are and respectively, which means that the Hopf bifurcation is supercritical and the periodic orbits are born stably, that is to say the system exists a stable limit cycle when (or ). And we can see that when (), the system (3.5) only exists the disease-free equilibrium, which is globally asymptotically stable, shown in Fig.1(b); when () and (), the system (3.5) exists an unstable disease-free equilibrium and a globally stable endemic equilibrium, shown in Fig. 1c and (e); when (), the system (3.5) exists an unstable disease-free equilibrium and an unstable endemic equilibrium, but exists a bifurcated periodic solution around this unstable endemic equilibrium, shown in Fig. 1d.

Fig. 1

a: Bifurcation diagram in () plane with (associated with parameter ), where , , , , , , , , , . The blue curve represents stable equilibrium, the red dashed curves represent the existence of bifurcated periodic solution around the equilibrium and a black dot marked with “Hopf” ( () and ()) indicates that the system (3.5) goes through Hopf bifurcation. (b-e): The time series of the system (3.5) with different , and in (b-e), is 0.4950, 2.4752, 14.8511, 34.6526 (the corresponding values of are 0.1, 0.5, 3, 7) and the initial values are (1609, 80, 52, 90), (), (2901.6825, 36.9450, 61.3725, 73.8900), (6841.3621, 50.0298, 100.6080, 100.0596), respectively

A backward bifurcation is shown in Figs. 2, and 2a illustrates the saddle-node bifurcation at and backward bifurcations at . Fig. 2b indicates that the system (3.5) has two endemic equilibria for , denoted as , where and , since functions and have two intersections for in such scenario. When (), the system (3.5) exists a unique globally stable disease-free equilibrium, shown in Fig. 1c. As parameter increases the DFE coexists with a stable endemic state for (i.e., ), shown Fig. 1d, and as further increasing parameter leads to system (3.5) stabilizes to the globally stable endemic equilibrium for (), shown in Fig. 2e.

Fig. 2

a: Bifurcation diagram in () plane with , where , , , , , , , , . The blue (blue dashed) curves represent stable (unstable) equilibrium and a black dot marked with “SN” () indicates that the system (3.5) goes through saddle-node bifurcation. b: The curves of functions and with . It shows that functions and have two intersections for , so system (3.5) has two endemic equilibria for . (c-e): The time series of the system (3.5) with different parameters , and in c-e, the values of parameter are 0.5, 3, 5 respectively

The forward-backward bifurcation is illustrated in Fig. 3a where saddle-node bifurcation is denoted by “SN” as a black dot. It is interesting to note that three endemic equilibria (denoted as , where and ) are feasible, shown in Fig. 3c-d, where functions and have three intersections for . It follows that when (i.e., ), the system (3.5) only exists the disease-free equilibrium, which is globally stable; when (i.e., ) and (), the system (3.5) has a stable endemic equilibrium; when (i.e., ), the system (3.5) has an unstable endemic equilibrium and two stable endemic equilibria, shown in Fig. 3b. It follows from Fig. 3b that and are locally stable, then solutions of I(t) initiating from near the equilibrium converge to either or .

Fig. 3

a: Bifurcation diagram in () plane with , where , , , , , , , , , the blue (blue dashed) curves represent stable (unstable) equilibrium and a black dot marked with “SN” ((), ()) indicates the system (3.5) goes through saddle-node bifurcation. (b): The time series of the system (3.5) for with different initial value. (c): The curves of functions and . (d): A partial enlargement of subplot (c)

Another pattern of the forward-backward bifurcation is illustrated in Fig. 4a where saddle-node bifurcation occurs for and Hopf bifurcation may also occur, and Fig. 4b is a partial enlargement of subplot Fig. 4a. We find two Hopf bifurcation points at and and two saddle-node bifurcation points at and . The first Lyapunov coefficients at those two Hopf points are both positive, which means that the Hopf bifurcations are subcritical and the bifurcated periodic orbits are born unstable, shown in Fig. 5, where some orbits (blue curves) initiating from near the equilibrium approach the stable endemic equilibrium (represented by red dot), while others (black curves) starting the value far away the equilibrium leave away the equilibrium, and consequently there is an unstable limit cycle between two kinds of orbits. To further analyze how periodic orbits change, we plot local bifurcation diagrams (Fig. 4c and d) and the variation in period of cycles versus (Fig. 4e and f). Figure 4c or d shows that once the periodic orbit appears, the minimum value of I(t) of the periodic orbit is infinitely close to the value of I(t) at the unstable saddle point (or ) at and . Figure4e and f shows that the periods eventually go to infinity at and , at which the homoclinic bifurcations occur and the limit cycle is replaced by a homoclinic orbit [30].

Fig. 5

Local phase portraits of system (3.5) with different ( in a, in b), other parameters are the same in Fig. 4a. The black dots represent the initial points and red dots represent stable equilibria. Here some orbits (blue curves) initiating from near the equilibrium approach the stable endemic equilibrium, while others (black curves) starting the value far away the equilibrium leave away the equilibrium, and consequently there is an unstable limit cycle between two kinds of orbits

In summary, we have illustrated various bifurcations including forward bifurcation, backward bifurcation, and two patterns of forward-backward bifurcations with saddle-node bifurcation and Hopf bifurcations. Theoretical analysis indicates that three or four positive equilibria are feasible for , but unfortunately we can not give the numerical illustrations.

In order to analyze the influence of media coverage on occurrence of backward bifurcation, we setwhich can be regarded as a parameter to measure the effect of psychological impact of media reported numbers of infectious individuals. On the one hand, considering the expression of function , we have . Algebraic calculation yieldswhich implies that we can increase k to make be less than 0, and hence multiple equilibria may be infeasible for . On the other hand, because of , we can increase k to make hold true to avoid the occurrence of backward bifurcation when on the basis of Theorem 1. Recall that when and , the system (2.1) may occur backward bifurcation and the basic reproduction number no longer act as the threshold to distinguish whether the disease dies out or not, but we can improve the effectiveness of psychological impact of media reported numbers of infectious individuals to avoid the occurrence of backward bifurcation, and consequently eliminating the infectious diseases becomes relatively easy.

In general, when only considering nonlinear recovery function, system (3.5) may have multiple endemic states due to backward bifurcation and forward-backward bifurcations. In particular, the system can exhibit rich dynamic behaviors including the existence of multiple endemic equilibria, Hopf bifurcation, saddle-node bifurcation, homoclinic bifurcation and unstable cycle. Comparing with the corresponding system with linear recovery function [4, 8, 12] we can conclude that nonlinear recovery may induce multiple stable endemic states, coexistence of the stable periodic solution (or endemic state) and the disease-free equilibrium.

We then investigate whether nonlinear media response induces those complex dynamics in the next section by briefly analyzing the case when considering system (2.1) with only saturated media coverage.

Considering only saturated media coverage in system (2.1), we have the following model equationsThe existence and stability of disease-free equilibrium for system (4.18) are same as those for the system (2.1). We let be the possible endemic equilibrium, then it satisfiesDenoteAlgebraic calculation yieldsfor . Further, we can testify thatif and only if andwhere . Combining the monotonicity of , we have function has no zero solution on the interval when , and there exists a unique I such that , denoted by (), for . However, the concrete form of cannot be represented.

According to the above analysis, we may come to the following conclusion.

Theorem 7

When , the system (4.18) has a unique endemic equilibrium, denoted by , which satisfies (4.19) and ; when , the system (4.18) has no endemic equilibrium.

Similarly we can investigate the global stability of the disease-free equilibrium and this unique endemic equilibrium when it exists. In the following we only give the main conclusions and the detailed proof processes are given in Appendix B.

Theorem 8

The disease-free equilibrium of the system (4.18) is globally asymptotically stable for . The unique endemic equilibrium of the system (4.18) is locally asymptotically stable for , and it is globally asymptotically stable if , and .

When only considering the saturated media coverage, we find that the basic reproduction number actually acts as the threshold to distinguish whether the disease dies out or not. Although we get the global stability of the unique endemic state under and certain conditions, numerical simulations imply that the system (4.18) has no other complex dynamics for . These results are also consistent with the conclusion in [12], indicating that saturated media coverage may not induce more complicated dynamics, compared with the corresponding model with linear media coverage. However, it is easy to testify that for , implying that saturated media impact does affect the disease prevalence. Therefore, the saturation phenomenon of media coverage may hardly influence on the dynamics in terms of inducing multiple endemic equilibria or bifurcations, but does affect the endemic infection level.

A case study and sensitivity analysis

In this subsection we initially parameterize the proposed model with surveillant data and data related to news items during early stage of COVID-19 outbreak in mainland China. We obtained the reported cumulative number of confirmed COVID-19 cases in China from the National Health Commission of the People’s Republic of China [31]. Although the first case was reported in December 2019, a new confirmed case was not reported until 10 January 2020, we used data from 10 January to 29 January, 2020, as shown in Fig 6a. Note that the case data was released and analyzed anonymously. We also obtained daily weighted average number of media items from 7 major websites during January 10-29, 2020, as in reference [11], which is shown in Fig. 6a.

Fig. 6

a :The reported cumulative number of confirmed cases and the average daily number of media items from January 10 to 29 . b and c : Data fitting for the data from January 15 to 29, 2020. The circles in b and c represent the cumulative number of confirmed cases, the average daily number of media items, from January 10 to 29 respectively. The curves are the best fitting curves of model (2.1) to these data

During the initial stage of the outbreak, we assume that the total population to be Wuhan residents, the number of infected population was the reported confirmed population and no individual was recovered. We use the Least Square Method to fit model (2.1) to data and the fitting results are shown in Fig. 6b and c, and the estimated parameter values are listed in Table 2. Based on the above-mentioned parameter estimations, we then calculate the basic reproduction number as 2.8611. Note that this estimation is within the confidence interval based on likelihood-based methods of [32, 33], while it is less than those estimated by the dynamic models without considering media impact [34, 35]. This indicates that media coverage may lead to the new infections decline by influencing individuals’s behaviour changes. On the basis of the estimated parameter values, we show the dynamics of proposed model in Fig. 7. It shows that the disease eventually converges to an equilibrium with a great number of infected individuals. Note that this seriousness of disease infection is because we do not consider other prevent and control interventions.

Table 2

Estimated initial values of variables and parameters for system (2.1)

Variables	Description	Initial value	Resource
S(0)	Initial susceptible population	1.1081×10^7	LS
E(0)	Initial exposed population	72.6413	LS
I(0)	Initial infected population	41	Data
M(0)	Initial Media items	16.3	Data

Fig. 7

The real dynamics of proposed model on the estimated real parameter values, where b is a partial enlargement of a

In order to quantify the influence of different parameters on the uncertainty in the model output, especially for the cumulative number of the infected individuals on February 10, 2020 ( i.e., at the end of the first month of the outbreak), we further conducted global sensitivity analysis on model (2.1). Basically, there are two main sensitivity analysis techniques: local and global methods. Global sensitivity analysis (GSA) techniques are more appropriate if parameter values have large uncertainties or models are nonlinear [36]. Here, two common GSA methods, the Morris method and the method of Sobol, are carried out. For Morris method [37, 38], there are two sensitivity measures that can be calculated for each input factor. One is the mean elementary effect (), which estimates the overall effect of the parameters on the output, and the other is the standard deviation of the elementary effects (), which estimates the combined effects of a parameter, including nonlinearity and effects due to interactions with other parameters. For the method of Sobol [37], there are also two sensitivity measures that can be calculated. One is the first-order sensitivity index (), which represents the main effect (direct contribution) of each parameter to the variance of the output, and the other is the total effect index (), which accounts for the total contribution to the output variation due to a parameter [37]. Note that the absolute values of the elementary effect () is further derived [39], and has been proved to be a good proxy for the total effect index, , based on variance-based measures [40].

For each method, 10 parameters were allowed to vary in a uniform range of +50% of their estimate values (shown in Table 2) and the effect of population size in initial epidemiological state (S(0), E(0), I(0) and M(0)) was not considered. We focused on the effect of each parameter on the cumulative number of the infected individuals at the end of the first month of the outbreak. The sensitivity coefficients obtained by the method of Sobol are given in Table 3. The sum of first-order sensitivity indices is approximately 0.6501 and the sum of total effect indices is 1.4150. The indicates that there should be interaction among model parameters. The first three parameters that have great influence on the output of the model are , and , which are all related to the characteristics of the disease itself. Note that parameter is also related to contact rate. Parameters , , and have the next greatest influence on the output results (their influence decreases in order). The first-order sensitivity index and total effect index for parameters and d are the smallest, which means that they have the least impact on the cumulative number of the infected individuals. It makes sense since the demographic process, due to its long timescale in comparison to short epidemic timescale, hardly influences on the transmission dynamics. In addition, total effect indices of parameters other than , d and are much greater than their first-order effect indices, indicating that these parameters interact with other parameters.

Table 3

The sensitivity coefficients obtained by the method of Sobol

Parameters	Λ	β	α	d	γ
The first-order sensitivity index	0.0024	0.3800	0.0097	0.0024	0.0147
The total effect index	0.0002	0.7218	0.0297	0.0002	0.0371

The similar conclusion can be obtained by Morris method, shown in Table 4. There is a slight difference, that is, the effect of parameter is much greater than that of parameter d and . Parameters other than , d and have high levels of , indicating there must be interaction among model parameters as well [37].

Table 4

The sensitivity coefficients obtained by Morris method

Parameters	Λ	β	α	d	γ
μ^∗	0.0755	196246.5489	22794.3909	0.1061	32749.9913
σ	0.0894	221419.2251	31416.1966	0.1278	26939.2488

It is worth noting that many studies have demonstrated that media coverage can be considered an effective way to mitigate the spread of disease during the early stages of an outbreak [8, 11, 19, 20]. To further examine the impact of saturation of media coverage and recovery on disease infections, we plotted the cumulative number of the infected individuals during the initial month of the outbreak with various values of , , and , as shown in Fig. 8. It shows that the less saturation of the media coverage () and recovery (), or the larger media reporting rate () or weight of media effect on contact rate sensitive to media items () is, the lower the cumulative number of the infected individuals during the initial month of the outbreak. It indicates that providing adequate medical resources and improving media response to infection and individuals’ response to mass media may benefit disease control during the early stage of the outbreak.

Fig. 8

The cumulative number of infected persons during the initial month of the outbreak with different values of , , and , respectively

Conclusion and discussion

During the early stage of emerging and re-emerging infectious diseases (say, COVID-19 infection) massive news coverage significantly influences the disease spread by generating psychological/behavioural changes, while shortage of media recourses worsens the disease infection. In this study, we proposed a mathematical model with an extra compartment (M(t)) to represent the level of media coverage, and further include double nonlinear functions to investigate the impacts of saturated media coverage and limited medical resources on disease transmission. The model provides a natural description for real transmission of an emerging infectious disease given in practical situations. Given great difficulties caused by high nonlinearity in theoretical analysis, we separately considered subsystems with only nonlinear recovery (i.e., ) or with only saturated media impact (i.e., ). For a special case, we technically constructed a new Lyapunov function using compound matrices and geometric ideas [4, 8, 12] to prove the global stability of the endemic state for the high-dimension model. Note that this newly constructed Lyapunov function can be generalized to other models.

We initially investigated the subsystems with only nonlinear recovery (i.e., ). Theoretical analysis shows that backward bifurcation occurs under certain parameters, and multiple equilibria can coexist in such system. And numerical simulations reveal some rich dynamic behaviors, including forward-backward bifurcation, Hopf bifurcation, saddle-node bifurcation, Homoclinic bifurcation and unstable limit cycle. We further consider the other special case where the recovery term is linear (i.e., ) to clearly investigate the contribution of only saturated media coverage on those complex dynamics, we proved the subsystem (4.18) has a unique endemic equilibrium which is locally stable for and a stable disease-free equilibrium for otherwise. The global stability of the endemic state under certain conditions is also obtained.

We noted that the subsystem (4.18) with saturated media impact exhibits the similar dynamics to the model with linear media impact, and further the basic reproduction number actually acts as the threshold to distinguish whether the disease dies out or not. We then conclude the saturated media impact slightly affects the dynamics. While including media coverage as an independent variable, nonlinear recovery may induce backward bifurcation and other rich dynamics such as Hopf bifurcation, Homoclinic bifurcation and unstable limit cycle. This means a combination of nonlinear recovery and media coverage significantly affect disease transmission in terms of causing oscillations. We also numerically analyzed effect of key parameters on occurrence of backward bifurcation, and obtained that increasing the effect of psychological impact of media reported numbers of infectious individuals can make the backward bifurcation less likely. Then improving media response to infection and individuals’ response to mass media may be beneficial for disease control.

Further, we parameterized the proposed model on the basis of the COVID-19 case data and data on news items in mainland China, and estimated the basic reproduction number to be 2.86. Sensitivity analyses were carried out to quantify the relative importance of parameters in determining the cumulative number of the infected individuals during the early stage of the outbreak. Combining with numerical analyses, we can conclude that providing adequate medical resources and improving media response to infection or individuals’ response to mass media may benefit disease control during the beginning of the COVID-19 outbreak. The model still has some shortcomings, such as not considering the immunity to disease of recovered individuals, which will be the content of further work.