key: cord-0703049-1oeht7fi
authors: Li, Tangjuan; Xiao, Yanni
title: Complex dynamics of an epidemic model with saturated media coverage and recovery
date: 2022-01-16
journal: Nonlinear Dyn
DOI: 10.1007/s11071-021-07096-6
sha: 8a362b7773595dc1c7c5e42c73b2745d2837752c
doc_id: 703049
cord_uid: 1oeht7fi

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.

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] [3] [4] [5] [6] . Some proposed the on-off media functions [5] [6] [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] [9] [10] [11] [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.

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 σ d , and then recover with a saturated function γ I 1+h 1 I and back into susceptibles class, where γ is the recovery rate of the infective individuals and h 1 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 (σ d E) with a saturated function δ 1+h 2 σ d E , where δ represents the reporting rate and h 2 is the parameter that measures the saturation effect of media reports. We further assume that media reports wane at the rate of μ d . 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 e −α M(t) , to describe the reduction factor in the incidence. Thus the model equations are

where is the rate of flowing into the population and d is the natural death rate. All parameters are nonnegative constants. We can testify that the feasible region of model (2.1) is

The disease-free equilibrium E 0 = ( d , 0, 0, 0) 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 is

where R 0 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.

Proof The Jacobian matrix concerned with the linearization of system (2.1) at E 0 is

So, the corresponding characteristic equation is

This implies that the disease-free equilibrium E 0 is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1. When R 0 = 1, zero is a simple eigenvalue of J E 0 and all other eigenvalues of J E 0 have negative real parts. Hence E 0 is a nonhyperbolic equilibrium, and the linearization cannot determine its stability. Here, we use center manifold [22] to analyze its stability. Let

We denoted the left and right eigenvalues of J E 0 as v and w respectively,

where m > 0. One can further choose m and n such that mn > 0 and vw = 1. Letγ = β d σ d d+σ d − d and φ =γ − γ , then R 0 < 1 if and only if φ < 0. According to Theorem 4.1 in [23] , the local dynamic behavior near R 0 = 1 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 f 2 and f 3 at E 0 when R 0 = 1 and φ = 0. Algebraic calculations show that

All other partial derivative of f 2 and f 3 are zero. So, we can get that

σ d d must hold, then we can get a < 0 if and only if

So we get the bifurcation at R 0 = 1 is backward when h 1 > h 1 ; the bifurcation at R 0 = 1 is forward when h 1 > h 1 , that is to say when R 0 = 1, E 0 is locally asymptotically stable if h 1 < h 1 and unstable if h 1 > h 1 . In particular, when h 1 = 0, h 1 < h 1 holds true, so the bifurcation at R 0 = 1 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 R 0 = 1, which is mainly associated with parameter h 1 , instead of parameter h 2 . 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 R 0 > 1 by using the standard methods [24] , shown in the following theorem (see Appendix A.1 for details).

Theorem 2 If R 0 > 1, there exists a positive constant η 0 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 R 0 = 1 is mainly associated with parameter h 1 , instead of parameter h 2 , we initially consider a special case: model with only nonlinear recovery (i.e., h 2 = 0), and conclude the dynamics of the proposed system theoretically and numerically.

When only considering the saturated recovery (i.e., h 2 = 0), the corresponding model equations are

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

3.1 Existence of the endemic equilibrium

As for system (3.5) , the possible endemic equilibrium E * * = (S * , E * , I * , M * ) satisfies

Denote A(I ) = d + γ 1+h 1 I , and

,

Then the existence of the solution of the last equation in (3.6) can be translated into the existence of the intersection of the functions f 1 (I ) and f 2 (I ). Algebraic calculation yields

Further, we can testify that

if and only if R 0 > 1 and

When (d +γ )(d +γ +σ d ) > σ d h 1 γ d , that is to say

we have B(I ) > 0 for ∀I > 0, so f 1 (I ) < 0 must holds for ∀I > 0 in this case, which implies that function f 1 (I ) decreases monotonically on [0, ∞) for h 1 <ĥ 1 .

Considering the zero point theorem and combining the monotonicity of f 2 (I ) for ∀I > 0, we conclude that function f 1 (I ) and f 2 (I ) have a unique intersection on the interval (0, d ) for h 1 <ĥ 1 and R 0 > 1; function f 1 (I ) and f 1 (I ) have no any intersection on the interval (0, d ) for h 1 <ĥ 1 and R 0 < 1.

, we get f (Î ) < 0 and f (I ) < 0 for I >Î , then f (I ) < 0 must hold true for I ≥Î , so we just need to consider if functions f 1 (I ) and f 2 (I ) have an intersection on [0,Î ]. For R 0 > 1, f (0) > 0 must hold, which implies that the system (3.5) at least has one endemic equilibrium for R 0 > 1 on the basis of the zero point theorem. i.e., the system has at least two positive equilibriums.

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

If h 1 <ĥ 1 , the system (3.5) has a unique endemic equilibrium for R 0 > 1, denoted by E 1 = (S * , E * , I * , M * ), which satisfies (3.6) and 0 < I * < I ; but it has no endemic equilibrium for R 0 < 1. If h 1 > h 1 , the system (3.5) exists at least one endemic equilibrium for R 0 > 1; the system (3.5) has no endemic equilibrium for R 0 < 1 and f max < 0, but it at least exists two endemic equilibrium for R 0 < 1 and f max > 0.

When h 1 <ĥ 1 and R 0 > 1 , the unique endemic equilibrium of system (3.5) is locally asymptotically stable.

Proof The Jacobian matrix concerned the linearization of system (3.5) at E 1 is

Let r (I * ) = γ 1+h 1 I * , we can get

Then the characteristic equation yields

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

When h 1 <ĥ 1 , the inequality

holds true sine the right side of inequality (3.13) is less than 0 if h 1 <ĥ 1 and the left side of inequality (3.13) must be greater than 0. And the inequality (3.13) is equivalent to the inequality

where r = γ 1+h 1 I * . The inequality (3.14) is further equivalent to the following inequality

which is naturally satisfied. Then in this case, a i > 0 (i = 2, 3) and a 1 a 2 − a 3 > 0 must hold true from the expressions of a i (i = 2, 3) and a 1 a 2 − a 3 . So if h 1 <ĥ 1 , a i > 0 (i = 1, 2, 3) and a 1 a 2 −a 3 > 0 hold true. According to the Hurwitz criterion, the roots of the polynomial λ 3 +a 1 λ 2 +a 2 λ+a 3 = 0 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 h 1 <ĥ 1 . That completes the proof.

The disease-free equilibrium E 0 of system (3.5) is globally asymptotically stable for h 1 <ĥ 1 and

Proof Let's consider a Lyapunov function

For h 1 <ĥ 1 , we can obtain 

That completes the proof.

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

Proof The Jacobian matrix of system (3.5) reads

The third additive compound matrix of J is

and the corresponding linear compound system reads

where

Then 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 as

where ι 1 ,ι 2 will be determined later. It is easy to show the system (3.5) admits a unique nonnegative bounded solution for all t ≥ 0 and any non-negative initial value.

Combing the Theorem 2, there exist two positive constants C 1 , C 2 such that

Recalling thaṫ

Then, discussing the derivative of V along the trajectory of the linear compound system (3.16) for different conditions and denoting D + the right-hand (total) derivative with respect to t.

Case (1):

Case (2):

Case (3):

Case (4): [26] , E 1 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 h 1 >ĥ 1 on the basis of Theorem 3. Here we have only analyzed a few specific cases, and others are similar but complicated.

Let's initially consider the case when h 1 >ĥ 1 and R 0 > 1. Algebraic calculation yields the quadratic equation B(I ) = 0 has a positive real root, denoted by I 1 , and a negative real root,

Then we have f 1 (I ) > 0 for 0 < I < I 1 and f 1 (I ) < 0 for I 1 < I < d . Further, we have 3 < 0, we then obtain that the function f 1 (I ) is a strictly convex function that increases on the interval (0, I 1 ] and decreases on the interval (I 1 , d ). Similarly we can analyze the concavity of the function f 2 (I ), which is given Case F : (A.6) holds and Q > 0

1 or 2 or 3 (See Table  6 for details)) 0 or 1 or 2 or 3 or 4

.. and H are given for discussing concavity of function f 2 (I ), referring to Table 5 in Table 5 by considering eight cases (case A-H), see Appendix A.2 for details.

Combining the monotonicity and concavity of functions f 1 (I ) and f 2 (I ), we discuss the number of local equilibria in the following three scenarios according to the number of intersection of functions f 1 (I ) and f 2 (I ) in the interval of (0, d ] (Table 1) . Scenario 1): If cases B, E or F hold, that is to say, f 2 (I ) is a strictly concave function of monotonically increasing. And the function f 1 (I ) is a strictly convex function that increases on the interval (0, I 1 ] and decreases on the interval (I 1 , d ). Combining (3.7) and (3.8), we can conclude the system (3.5) has a unique endemic equilibrium for R 0 > 1. When case A holds, f 2 (I ) 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

(I * 31 < I * 32 < I * 33 ) respectively. As for the case D, we can do the same analysis.

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

As for the case when h 1 >ĥ 1 , R 0 < 1, 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.

A forward bifurcation is illustrated in Fig.1 , and Fig.1(a) gives variation in equilibrium level of infections (I * ) with the basic reproduction number (associated with parameter ) and corresponding bifurcations, say, R 0 = 5.0738 ( = 1.0249) and when R 0 = 5.0738 ( = 6.0878), 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 −7.9051 × 10 7 and −2.2751 × 10 7 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 5.0738 < R 0 < 30.1368 (or 1.0249 < < 6.0878). And we can see that when R 0 < 1 (0 < < 0.2020), the system (3.5) only exists the disease-free equilibrium, which is globally asymptotically stable, shown in Fig.1(b) ; when 1 < R 0 < 5.0738 (0.2020 < < 1.0249) and R 0 > 30.1368 ( > 6.0878), the system (3.5) exists an unstable disease-free equilibrium and a globally stable endemic equilibrium, shown in Fig. 1c and (e); when 5.0738 < R 0 < 30.1368 (1.0249 < < 6.0878), 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 . A backward bifurcation is shown in Figs. 2, and 2a illustrates the saddle-node bifurcation at R 0 = 0.1532 and backward bifurcations at R 0 = 1. Fig. 2b indicates that the system (3.5) has two endemic equilibria for

, since functions f 1 (I ) and f 2 (I ) have two intersections for I > 0 in such scenario. When R 0 < 0.1532 ( < 0.7115), the system (3.5) exists a unique globally stable disease-free equilibrium, shown in Fig. 1c . As parameter increases the DFE E 0 coexists with a stable endemic state E 22 for 0.1532 < R 0 < 1 (i.e., 0.7115 < < 4.6453), shown Fig. 1d , and as further increasing parameter leads to system 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 E 3i = (S * 3i , E * 3i , I * 3i , M * 3i ), where i = 1, 2, 3 and I * 31 < I * 32 < I * 33 ) are feasible, shown in Fig. 3c -d, where functions f 1 (I ) and f 2 (I ) have three intersections for I > 0. It follows that when R 0 < 1 (i.e., 0 < < 2.5909), the system (3.5) only exists the disease-free equilibrium, which is globally stable; when 1 < R 0 < 1.0218 (i.e., 2.5909 < < 2.6474) and R 0 > 1.6565 ( > 4.2918), the system (3.5) has a stable endemic equilibrium; when 1.0218 < R 0 < 1.6565 (i.e., 1.0218 < < 4.2918), the system (3.5) has an unstable endemic equilibrium and two stable endemic equilibria, shown in Fig. 3b . It follows from Fig. 3b that E 31 and E 33 are locally stable, then solutions of I (t) initiating from near the equilibrium E 32 converge to either I * 31 or I * 33 . Another pattern of the forward-backward bifurcation is illustrated in Fig. 4a where saddle-node bifurcation occurs for R 0 < 1 and Hopf bifurcation may also occur, and Fig. 4b is a partial enlargement of subplot Fig. 4a . We find two Hopf bifurcation points at R 0 = 0.782390 and R 0 = 1.001113 and two saddle-node bifurcation points at R 0 = 0.752654 and R 0 = 1.417937. The first Lyapunov coefficients at those two Hopf points are both positive, which means that the Hopf bifurcations are subcritical and the bifur- 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 sta-ble 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 cated 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 R 0 ( 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 E * 21 (or E 0 ) at R 0 = 0.783296 and R 0 = 1.001112. Figure4e and f shows that the periods eventually go to infinity at R 0 = 0.783296 and R 0 = 1.001112, at which the homoclinic bifurcations occur and the limit cycle is replaced by a homoclinic orbit [30] .

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 R 0 < 1, but unfortunately we can not give the numerical illustrations.

In order to analyze the influence of media coverage on occurrence of backward bifurcation, we set

which 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 f 2 (I ), we have f 2 (I ) = d d+σ d σ d β e k A(I )I . Algebraic calculation yields

which implies that we can increase k to make f max be less than 0, and hence multiple equilibria may be infeasible for R 0 < 1. On the other hand, because of ∂ h 1 ∂k > 0, we can increase k to make h 1 < h 1 hold true to avoid the occurrence of backward bifurcation when R 0 = 1 on the basis of Theorem 1. Recall that when h 1 >ĥ 1 and f max > 0, 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 occur-rence 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 forwardbackward 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 equations

(4.18)

The existence and stability of disease-free equilibrium for system (4.18) are same as those for the system (2.1).

We let E 1 = (S * , E * , I * , M * ) be the possible endemic equilibrium, then it satisfies Denote

Algebraic calculation yields

for ∀I > 0. Further, we can testify that

if and only if R 0 > 1 and

Combining the monotonicity of η(I ), we have function η(I ) has no zero solution on the interval [0, ∞) when R 0 ≤ 1, and there exists a unique I such that η(I ) = 0, denoted by I * (0 < I * < I ), for R 0 > 1. However, the concrete form of I * cannot be represented.

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

Theorem 7 When R 0 > 1, the system (4.18) has a unique endemic equilibrium, denoted by E 1 = (S * , E * , I * , M * ), which satisfies (4.19) and 0 < I * < I ; when R 0 ≤ 1, the system (4.18) has no endemic equilibrium.

Similarly we can investigate the global stability of the disease-free equilibrium E 0 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. 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 R 0 > 1 and certain conditions, numerical simulations imply that the system (4.18) has no other complex dynamics for R 0 > 1. 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 ∂ I * ∂h 2 > 0 for R 0 > 1, 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.

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 .

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 R 0 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.

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 (S i ), which represents the main effect (direct contribution) of each parameter to the variance of the output, and the other is the total effect index (S T i ), 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, S T i , 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- 

Average number of daily media items Fitting curve 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 β, σ d and γ , which are all related to the characteristics of the disease itself. Note that parameter β is also related to contact rate. Param-eters α, δ, h 2 and h 1 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 com- Table 4 . There is a slight difference, that is, the effect of parameter μ d is much greater than that of parameter d and . Parameters other than , d and μ d have high levels of σ , indicating there must be interaction among model parameters as well [37] .

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 α, δ, h 1 and h 2 , as shown in Fig. 8 . It shows that the less saturation of the media coverage (h 2 ) and recovery (h 1 ), 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.

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., h 2 = 0) or with only saturated media impact (i.e., h 1 = 0). 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., h 2 = 0). Theoretical analysis shows that backward bifurcation occurs under cer-tain 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., h 1 = 0) 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 R 0 > 1 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.

Funding This research was funded by the National Natural Science Foundation of China grant numbers: 12031010, 11631012(YX).

Availability of data and material Not applicable. 

Then X is the positive invariant set of X 0 , and X 0 is relatively closed in X . In order to prove that the disease is uniformly persistent, it suffices to show that ∂ X repels uniformly the solutions of (2.1) in X 0 .

First, we can see the system (2.1) is point dissipative since the attracting domain = {(S, E, I, M) ∈ R 4 + : Since

Then we can prove that for any solution u(t, x 0 ),

is the stable manifold of E 0 . Let's prove it by proof by contraction. If lim sup t→∞ u(t, x 0 ) − E 0 < δ 0 , i.e., there is a positive constant η 0 , for any t > η 0 , such that

then when t > η 0 , we have

The characteristic equation of the Jacobian matrix of (A.1) at (0, 0) is

which implies that the characteristic equation has a positive root, then any positive solution of (A.1) tends to infinity as t tends infinity, this contradicts 0 < S + E < I < d . Thus, E 0 is an isolated invariant set in set X , and W S (E 0 ) ∩ X 0 = ∅.

In conclusion, each of the forward orbital in M ∂ converges to E 0 , and E 0 is aperiodic in M ∂ . By Theorem4.6 of [24] or Theorem4.3 and Remark 4.3 of [41] , we can conclude that ∂ X repeals uniformly the solutions of (2.1) in X 0 .

Algebraic calculation yields

where m = 2γ μ d h 1 αδ , and we only need to consider the case C > 1. We can get lim C→+∞ P 1 (C) − P 2 (C) = +∞

= Q Algebraic calculation yields P 1 (C) > 0, P 1 (C) > 0, P 2 (C) > 0 and P 2 (C) < 0, which means that P 1 (C) is a strictly concave function of monotonically increasing, and P 2 (C) is a strictly convex function of monotonically increasing for C > 1.

If Q > 0, we can show that functions P 1 (C) and P 2 (C) intersect at one point if and only if there is a constantĈ > 1 such that P 1 (Ĉ) = P 2 (Ĉ) and P 1 (Ĉ) = P 2 (Ĉ), that is equations

have a unique solution forĈ > 1. We can further testify that equations (A. Similarly, we can show that functions P 1 (C) and P 2 (C) intersect at two points if and only if there is a constant C > 1 such that P 1 (Ĉ) = P 2 (Ĉ) and P 1 (Ĉ) > P 2 (Ĉ), that is there is a constantĈ > 1 such that equations Let's call these two points of intersectionĈ 1 andĈ 2 (Ĉ 1 <Ĉ 2 ) in this case. We can further testify that when there is no intersection between P 1 (C) and P 2 (C).

Denoting C = 1 + h 1 I > 1, and considering the positive linear correlation between C and I and the expression for the function f 2 (I ), we can tell whether Table 5 The concavity of f 2 (I )

Conditions Conclusion If there is no intersection between P 1 (C) and P 2 (C), f 2 (I ) > 0 f or I > 0.

On the other hand, if Q < 0, combing the monotonicity and concavity of functions P 1 (C) and P 2 (C), we can testify functions P 1 (C) and P 2 (C) have a unique intersection for C > 1, denoted asĈ. Similarly, let I =Ĉ −1 h 1 , we can get Combing that 0 < I < d , we can get to the concavity of f 2 (I ), which is show in Table 5 or Table 1 . The number of endemic equilibrium of the case C is given in Table 6 .

Proof We first testify the global stability of the diseasefree equilibrium E 0 for R 0 ≤ 1, let's consider a Lyapunov function 

I * [1 + h 2 (d + r(I * ))I * ] 2 > 0,

It's easy to show that the characteristic equation (B.1) has a negative root −d, and other roots are the roots of the polynomial λ 3 + n 2 λ 2 + n 1 λ + n 0 = 0. Further, we have n 2 n 1 − n 0 = μ 2 d [g(I * ) + 2d + σ d + γ ] + μ d [2d + γ + σ d + g(I * )] 2 + [2d + γ + σ d + g(I * )](d + γ + σ d )g(I * )

According to the Hurwitz criterion, the roots of the polynomial λ 3 + n 2 λ 2 + n 1 λ + n 0 = 0 all have negative real parts, then the unique endemic equilibrium of the system (4.18) is locally asymptotically stable if it is feasible (R 0 > 1).

Finally, we use compound matrices and geometric ideas to develop global results under mild restrictions as well. The Jacobian matrix of system (4.18) reads 

Then 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 as

Modeling the effect of literacy and social media advertisements on the dynamics of infectious diseases

Analysis of a model of the glucose-insulin regulatory system with two delays

The effect of constant and pulse vaccination on sis epidemic models incorporating media coverage. Commun

The effects of media reports on disease spread and important public health measurements

A filippov system describing media effects on the spread of infectious diseases

Media impact switching surface during an infectious disease outbreak

Dynamics of an infectious diseases with media/psychology induced non-smooth incidence

Modeling and analysis of effects of awareness programs by media on the spread of infectious diseases

Analysis of an epidemic system with two response delays in media impact function

Media coverage and hospital notifications: correlation analysis and optimal media impact duration to manage a pandemic

Effects of media reporting on mitigating spread of covid-19 in the early phase of the outbreak

Optimal media reporting intensity on mitigating spread of an emerging infectious disease

Dynamics of infectious diseases: local versus global awareness

Saturation recovery leads to multiple endemic equilibria and backward bifurcation

Analysis of stability and bifurcation for an seir epidemic model with saturated recovery rate

Dynamic behaviors of a modified sir model in epidemic diseases using nonlinear incidence and recovery rates

Backward bifurcation and sensitivity analysis for bacterial meningitis transmission dynamics with a nonlinear recovery rate

Projecting hospital utilization during the covid-19 outbreaks in the united states

Effects of medical resource capacities and intensities of public mitigation measures on outcomes of covid-19 outbreaks

Media/psychological impact on multiple outbreaks of emerging infectious diseases

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Qualitative theory of second-order dynamic systems

Dynamical models of tuberculosis and their applications

Persistence under relaxed point-dissipativity (with application to an endemic model)

Dynamics of differential equations on invariant manifolds

Global stability of seirs models in epidemiology

A model of hiv/aids with staged progression and amelioration

Global dynamics of an seir epidemic model with vertical transmission

Matcont and cl matcont: Continuation toolboxes in matlab

An hiv infection model based on a vectored immunoprophylaxis experiment

Transmissibility of 2019-nCoV

Early transmission dynamics in wuhan, china, of novel coronavirus-infected pneumonia

Estimation of the transmission risk of the 2019-ncov and its implication for public health interventions

The reproductive number of covid-19 is higher compared to sars coronavirus

A methodology for global-sensitivity analysis of time-dependent outputs in systems biology modelling

Global Sensitivity Analysis: the Primer

Sensitivity analysis of infectious disease models: methods, advances and their application

An effective screening design for sensitivity analysis of large models

Sensitivity analysis of environmental models: a systematic review with practical workflow

Chain transitivity, attractivity, and strong repellors for semidynamical systems

Amount Range I 1 ≤Ĩ 1 <Ĩ 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 )OneTwoOneOne I * > I 1 rium E 0 of the system (4.18) is globally asymptotically stable if R 0 ≤ 1. That completes the proof. Then, we prove the stability of this unique endemic when it exists. According to the Theorem 7, the system (4.18) has a unique endemic equilibrium E 1 for R 0 > 1. The Jacobian matrix concerned the linearization of system (4.18) at E 1 isand the characteristic equation yields