Introduction

Emerging and reemerging infectious diseases including the 2003 severe acute respiratory syndrome (SARS) and the 2009 A/H1N1 influenza epidemic have become a major cause of mortality and morbidity in emergency situations. News reports have the potential to modify a community’s knowledge of emerging infectious diseases, and affect people’s attitudes and behaviours during infectious disease outbreaks [1, 2]. People informed by media reports can take precautions ranging from washing hands, wearing protective masks to avoiding social contact with infected individuals, to reduce their susceptibility. Informed infective individuals will also take measures to protect themselves from being exposed to others to reduce infectivity. It has been shown that behaviour change during infectious disease outbreaks can curb the effects of infectious diseases in populations.

In recent years, a growing number of studies have focused on understanding and quantifying the impact of such behaviour influencing factors on the spread of infectious diseases [4–18]. A number of studies have employed mathematical models to assess the impact of media reports on emerging infectious disease prevention and control [4, 5, 19–30]. Recently, Greenhalgh et al. presented a brief and nice commentary on the literature related to awareness and their effects on the dynamics of diseases. In summary, it has been found that there are three main methods being used to incorporate behaviour change in mathematical models due to awareness of disease. In the first method, the incidence rate of the disease is reduced by some factors that depend on the numbers of infected individuals, hospitalized individuals or exposed individuals, due to education about preventative knowledge of the disease through media coverage. The common choice of the reduction factors is a saturated [4, 6, 31, 32] or exponential [20, 24, 27] growth function. For example, Liu et al. incorporated an exponential decreasing factor β0 = βe−a1E−a2I−a3H into the transmission coefficient (with exposed (E), infectious (I), hospitalized (H)) to illustrate the possible mechanism for multiple outbreaks of SARS due to the psychological impact. Cui et al. used the general nonlinear incidence function μ1 − μ2

f(I) to represent the media and education impact on the spread of the infectious disease. In, the authors focused on simple endemic models by modelling the contact rate as a function of the available information on the present and the past disease prevalence.

In the second method, a separate compartment that effectively represents the level of awareness in the population is introduced, and individuals in the population can move from the unaware to aware compartments [7, 15, 16, 19, 21–23, 28]. For example, in, the authors proposed a mathematical model by inducing behavioural changes in the population through delineation of the susceptible class into unaware susceptible and aware susceptible subpopulations. and explicitly introduced distinct compartments for unaware and aware individuals in each of the disease states, and transitions between respective unaware and aware compartments took place at constant rates.

In the third method, a compartment representing the awareness program is incorporated [8, 18, 21, 25]. Yan and Tang et al. described the effects of media reports on population infection by modifying the transmission rate β following an exponential function βe−pM with M representing the level of media reports. Further, in, the media reporting is introduced as a separate compartment in a mathematical model and the susceptible population is divided into three awareness levels, each with a different infection rate. In, the authors considered the interaction of disease outbreak and media impact by formulating a susceptible-infected-hospitalized-recovered framework of population. By extension, susceptible and infected populations are subdivided into aware and unaware since individuals modified their behaviors to reduce their transmissibility and infectivity, and the dynamics of media reports was incorporated by considering how media was influenced by the numbers of infected and hospitalized individuals.

The majority of the mathematical modelling studies described above have incorporated media impact either in the disease transmission term or by dividing the susceptible population into subgroups with various awareness levels. However, the relationship between mass media and disease spread can be more complex than these models portray. On one hand, media reporting influences the public awareness of the disease and affects the effectiveness of prevention measures. On the other hand, the severity of the disease has an impact on the degree of mass media reporting. We’ve known that, in the first and third methods, media impact is modelled through the inclusion of a “media function”, which is proportional to the number of infected individuals and/or the level of media reports, to reduce the incidence rate through increased protective behaviour. However, it remains unclear as to whether awareness of the number of infections, or the awareness of media reports best modify individual behaviour during an infectious disease outbreak. This falls within the scope of this study.

Herewithin, we establish a mathematical model incorporating media reports as a separate compartment by considering how media is influenced by disease statistics (number of newly observed individuals). Disease progression is characterized by an SEIR model of which the transmission rate is modified by a media function affected by the media reports and also the number of infected individuals. The model can be recognized as the combination of the first method and the third method of modelling. In our model, we formulate the novel media function f(I, M, α1, α2) with αi (i = 1, 2) denoting the weight of infected individuals and media reports, and examine their effects on disease spread. Further, we investigate an optimal control problem in order to seek the optimal reporting intensity of information to minimize the number of infected individuals (and costs). We parameterize the proposed model on the basis of the 2009 A/H1N1 data in Shaanxi province of China, and estimate the basic reproduction number and other unknown parameters. A sensitivity analysis is conducted to identify model parameters that most affect the peak magnitude of the epidemic, as well as the total number of infections over the entire epidemic.

Model

We are interested in studying the effects of I and M on the outcomes of an infectious disease outbreak/epidemic. We therefore consider a mass media compartment M and a media effect function that depends on both the number of infected I and media reports M. Consider an SEIR model that incorporates a compartment of media programs M, in which the media impact on the human behaviour is reflected in the contact rate.

{dSdt=Λ−f(I,M,α1,α2)βSI−μS,dEdt=f(I,M,α1,α2)βSI−σE−μE,dIdt=σE−γI−μI,dRdt=γI−μR,dMdt=ρσE−δM.(1)

where S(t), E(t), I(t), R(t) represent the susceptible, exposed, infective, and recovered populations, respectively, M(t) represents the number of news items. Here, Λ is the birth rate, μ is the natural death rate, σ is the progression rate from the exposed to infective classes, and γ is the recovery rate. The propagation of information depends on the number of newly observed individuals (σE), and ρ represents the reporting rate of the newly observed individuals. It is assumed that δ represents the spontaneous disappearance rate of media. The baseline transmission rate without media effect is represented by β and f(M, I, α1, α2) is used to modify the transmission rate, which is induced by the media effect. Finally, α1, α2 are the weights of media effect sensitive to infectives and media items, respectively. All the parameters are non-negative.

It’s obvious that the media impact on the behavior of humans increases as the number of infected individuals increases or the intensity of media reports increases. Thus the term of media impact factor reflecting the behavior change f(I, M, α1, α2) is a decreasing function with respect to the number of infected individuals I and the media intensity M, it should satisfy the following assumptions:

∂f(I,M,α1,α2)∂I≤0forallI>0,∂f(I,M,α1,α2)∂M≤0forallM>0,f(0,0,α1,α2)=1,f(I,M,α1,α2)→0asI→∞orM→∞.(2)

Here, we choose two different media functions,

f1(I,M,α1,α2)=e−α1I−α2M,

and

f2(I,M,α1,α2)=11+α1I+α2M.

Optimal control

In system (1), the reporting rate is proportional to the newly reported number of infected individuals. It is natural to ask whether we should enhance the reporting rate to reduce the total infected number and minimize the cost of media reporting. Thus our main purpose is to minimize the total number of infective individuals as well as the cost required to reduce or enhance the media reporting intensity.

Consider the optimal control problem to minimize the objective functional

J(u)=∫t0T[AI(t)+B2u2(t)]dt(3)

subject to

{dSdt=Λ−βf(I,M)SI−μS,S(t0)=S0>0,dEdt=βf(I,M)SI−σE−μE,E(t0)=E0≥0,dIdt=σE−γI−μI,I(t0)=I0≥0,dRdt=γI−μR,R(t0)=R0≥0,dMdt=u(t)ρσE−δM,M(t0)=M0≥0.(4)

where the coefficients A and B/2 are positive. Here we assume that A = 1, and that B/2 is the weight associated with the control u(t). Note that u(t) is a Lebesgue measurable function on a finite interval [t0, tend], where 0 ≤ u(t) ≤ umax, umax > 1, and 0 ≤ u(t) < 1 represents reduction in reporting intensity, whereas 1 < u(t) ≤ umax represents enhancement in reporting intensity.

Parameter values

We use data from the 8th hospital of Xi’an in Shaanxi province to study the effects of media reports. The data are fully available in S2 File. The data include information on the daily number of hospital notifications in the 8th hospital from September 3 to 30, 2009. Parameter values for Eq (1) are informed by the literature and are further estimated through model fits to the hospital notification data, using the Least Square Method.

To ensure that our model estimates of the basic reproduction number, R0, are from the exponential growth phase of infection, we assume data from September 3-21. Also, so that we can compare the effects of I and M on the optimal control of media reporting, we assume a mass media compartment and a media function that depends both on I and M. Table 1 lists the best-fit parameters determined for model (1), without media impact f(I, M, α1, α2) = f0(I, M) = 1 and with two different media functions f1 = e−α1I−α2M and f2=11+α1I+α2M, respectively.

Note that, for completeness, we consider the extended dataset from September 3-30 in a sensitivity analysis. Also, we have provided model fits and parameter values when the mass media compartment M is not included in the model in S1 File Appendix C. However, as we are interested in comparing the effects of I and M in the current study, we consider the full model with the M compartment in the optimal control and sensitivity analysis.

Equilibrium

The basic reproduction number for system (1) can be calculated as

R0=βσΛμ(γ+μ)(σ+μ)(5)

easily using the next generation method or the survival function method (see for a review of this method and other methods that are commonly used). Note that the basic reproduction number is independent of the mass media compartment. Also note that it is not affected by the media function f(I, M). Therefore, it is the same as it would be in a model without media impact, which means that media coverage does not play a role in affecting this epidemic threshold. This has also been observed in previous studies [4, 6, 11].

Omitting the equation for R, the system (1) can be rewritten as a four dimensional model. Here, the disease free equilibrium is E0(Λμ,0,0,0), and the endemic equilibrium E*(S*,E*,I*,M*) should satisfy

{Λ−f(I*,M*)βS*I*−μS*=0,f(I*,M*)βS*I*−σE*−μE*=0,σE*−γI*−μI*=0,ρσE*−δM*=0,(6)

Simplifying gives

S*=Λμ−δ(σ+μ)ρσμM*,E*=δρσM*,I*=δρ(γ+μ)M*,(7)

and

f(I*,M*)(Λ−δ(σ+μ)ρσM*)=μ(σ+μ)(γ+μ)βσ.(8)

Since I* can be expressed by M*, we can rewrite f(I*, M*) as h(M*)=f(δρ(γ+μ)M*,M*). Denote h(M)=f(δρ(γ+μ)M,M), then h(M) is a decreasing function with respect to M. Let

g(M)=h(M)(Λ−δ(σ+μ)ρσM)−μ(σ+μ)(γ+μ)βσ.(9)

Then we have g(0)=Λ−μ(σ+μ)(γ+μ)βσ>0 when R0>1, g(ρσΛδ(σ+μ))=−μ(σ+μ)(γ+μ)βσ<0, and g′(M)=h′(M)(Λ−δ(σ+μ)ρσM)−δ(σ+μ)ρσh(M)<0 holds true for any M satisfying Λ−δ(σ+μ)ρσM>0. Thus there must exist one and only one positive root M*<ρσΛδ(σ+μ) that satisfies g(M*) = 0 if R0>1. Particularly, if we choose

f1(I,M)=e−α1I−α2M,

and

f2(I,M)=11+α1I+α2M,

we have

M*=ρσΛδ(σ+μ)−ρ(γ+μ)α1δ+α2ρ(γ+μ)LambertW(μ(α1δ+α2ρ(γ+μ))βδe(α1δ+α2ρ(γ+μ))σΛδ(σ+μ)(γ+μ)),

and

M*=ρμ(γ+μ)(R0−1)μ(α1δ+α2ρ(γ+μ))+βδ,

respectively. (Please find detailed definition of Lambert W function in paper [10, 37]).

Note that the disease free equilibrium E0(Λμ,0,0,0) is globally asymptotically stable if R0≤1, and unstable if R0>1. Meanwhile, the unique endemic equilibrium E*(S*,E*,I*,M*) exists if and only if R0>1 and it is locally asymptotically stable if it is feasible. For more information about the stability of the disease free equilibrium and the endemic equilibrium, see S1 File Appendix A.

Existence of optimal control

Denote the control set U={u(t):0≤u(t)≤umax,t0≤t≤tend,u(t)isLebesguemeasurable}. The existence of optimal control can be shown using the results from Theorem 4.1 in. We can easily verify the following properties:

The set of control and corresponding state variable is non-empty, which can be shown by the boundedness of solutions of system (4) using the results from Theorem 9.2 in.The control set U is closed and convex by definition.The right-hand side of the state system is bounded above by a linear function in the state and control, since the solutions are bounded, which determines the compactness needed for the existence of the optimal control.The integrand of the objective functional is convex on the control u(t), and there exists q1 > 0, q2 > 1 such that AI(t)+B/2u(t)2≥q1|u(t)|q2−q3, where we can choose q1 = B/2 and q2 = 2.

Then we have that, for the control problem (3) and (4), there exists an optimal control u*∈U such that min J(u) = J(u*) on the interval [t0, tend].

Theorem 1

There exists an optimal control u* that minimizes J(u) over

U. Moreover, there exists adjoint functions

{λS′(t)=βf(I,M)I(λS−λE)+μλE,λE′(t)=σ(λE−λI−ρu(t)λM)+μλE,λI′(t)=−A+βf(I,M)S(λS−λE)+βSI∂f(I,M)∂I(λS−λE)+γ(λI−λR)+μλI,λR′(t)=μλR,λM′(t)=βSI∂f(I,M)∂M(λS−λE)+δλM,(10)

with the transversality conditions

λS(tend)=0,λE(tend)=0,λI(tend)=0,λR(tend)=0,λM(tend)=0.(11)

The optimal control u* is given by

u*(t)=max{min{−λMρσEB,umax},0}.(12)

For detailed derivation of the optimal control for the control problem (3) and (4), see S1 File Appendix B.

Numerical simulation

During the 2009 A/H1N1 influenza pandemic, media coverage was used to spread precaution information about the disease, which influenced human behaviours [11, 25]. Using data extracted from the initial laboratory-confirmed cases of 2009 A/H1N1 that were admitted to the 8th hospital of Xi’an for the period 3rd September to 21st September, and the Least Square Method, we estimated the parameters shown in Table 1 and Fig 1, with R-square value being 0.9588, 0.9577, 0.9583, using three different media functions f0, f1, f2. Note that the goodness of fit is significant, but similar for each case. This is not unexpected as the data spans 19 days only. Also, however, as shown in our discussion of R0, during the early stages of an epidemic there is little effect due to mass media reports. Therefore it is reasonable that the model fits and parameter values are similar comparing the model with (f1, or f2) and without (f0) media.

Using the Akaike Information Criterion (AIC) for Least-Squares case, AIC=nlog(RSSn)+2k, where n is the sample size, k is the number of parameters, and RSS denotes residual sum of squares of fitted model, we obtain an AIC of 90.8918 for the model without media impact (i.e. the model with f0), 99.3676 for the model with f1, and 99.0981 for the model with f2. We note that the model without media impact has the lowest AIC and models that do not consider the mass media compartment M also fit the data well (see S1 File Appendix C). However, as we are interested in understanding the effects of mass media reports M and known infectives I to an individual in the population, we continue our study considering the model with the M compartment, and media functions f1 and f2. Considering the R-square and AIC, we conclude that model (1) with f2 fits the observed data better than the model with media function f1.

To show the sensitivity of the basic reproduction number R0 with respect to the time interval considered, we also estimated key epidemic parameters and initial conditions considering the time periods between September 3rd and 23rd, 25th, 28th, and 30th, respectively. Results are presented in Table 2 and shown in Fig 2. When different periods are considered, the reproduction number varies from 1.8715 to 2.0463. The results indicate that, for each period we consider, the model with both f1 and f2 can fit the observed data well. Note that we have chosen to consider the model fits using data from the 3rd to 21st of September in the analysis below. This is done to ensure that we estimate R0 during the exponential growth phase of the outbreak.

Using data from September 3-21, the basic reproduction number is estimated as R0 = 1.8248 without media impact (f(I, M, α1, α2) = f0 = 1), which is lower than the estimates for R0 = 1.8715 for media function f(I, M, α1, α2) = f1 and 1.8716 for media function f(I, M, α1, α2) = f2. These three values of R0 are all in agreement with the result in, in which the mean value of the basic reproduction number was estimated as 1.794 with 95% confidence interval [1.3858, 1.9091]. Again, we note that the similar R0 values reflect the fact that there is little impact of media reports in the early days of the epidemic.

We plot the epidemic curves of system (1) without media function (f0 = 1) and with media function f1, f2, in Fig 3, using the parameter values listed in Table 1. This figure shows that the number of infected individuals is greatly reduced when media impact is considered. For example, the peak magnitude is reduced from 2091 to 1225 (1189) when we use the media function f1 (f2), giving a reduction of 41.4% (43.1%). The total number of infected individuals over one year is also reduced, from 87974 to 74391 (73370). It also follows from Fig 3 that there is no obvious difference on the epidemic prevalence for system (1) with media function f1 or function f2, though the number of media items looks very different for the different media functions, which may be attributed to the small numbers of media items and the small value of the weight of media effects sensitive to the media reports. This is also confirmed in Fig 4(A) and 4(B) that media functions f1 and f2 are almost the same under their estimated parameters.

The contour plots of Fig 5 show the dependence of the peak magnitude and the total number of infections on the weight of media effects sensitive to infected individuals α1 and the weight of media effects sensitive to the media reports α2, using the two different media functions. With increases in α1 or α2, both the peak magnitude of the number of infected individuals and the total number of infections over a year decrease greatly, indicating that the media effect reduces outbreak severity. To identify key parameters that influence the disease infection dynamics, we use Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCCs) to examine the dependence of the peak magnitude and total number of infections on corresponding model parameters. It follows from Figs 6 and 7 that, despite the fact that the baseline value of the transmission rate β and the recovery rate γ are the most sensitive parameters to the peak magnitude and the total number of infections, parameters related to the media coverage α1, α2, ρ, δ can also significantly affect the results. In particular, increases in the weight of infective cases α1 and media reports α2 in the media function will significantly reduce the peak magnitude and total number infections. Also, decreases in the media reporting rate ρ, and increases in the media waning rate will lead to more severe outbreaks.

To access the effectiveness of enhancing the media reporting on the epidemic outbreak, we plot the variation in peak magnitude and total infections with the corresponding parameter ρ, as shown in Fig 8(A) and 8(B). It shows that increasing ρ by 4 times from baseline value (while keeping other parameters fixed) can reduce the peak magnitude from 1225 to 779 (decreased by 36.4%) for the media function f1 or 1189 to 687 (decrease by 42.2%) for the media function f2, and also can reduce the number of total infections from 74391 to 66783 for f1 or 73370 to 64102 for f2.

To further investigate what pattern of media report is optimal in minimizing the number of infected individuals and costs, we simulate the optimal control system (4) and obtain the optimal control. Here, we fix the parameter values as listed in Table 1 and employ the Forward-Backward Sweep method. It follows from Fig 9(A) that the optimal control is to continuously strengthen/increase the media reports at the beginning of an epidemic, when the disease starts to spread. A maximum level should then be maintained in times surrounding the peak number of infections, and then it should be slowly weakened as the infection reduces in the population. It is clear that the optimal control is similar for the two different media functions f1 and f2, however, the optimal media reporting intensity for system (4) with media function f2 is stronger than that for system (4) with media function f1 (i.e., u2*>u1*).

Figs 9(B) and 10 show the optimal epidemic curves under the optimal reporting intensity. These figures indicate that the optimal control significantly reduces the peak magnitude and total number of infected individuals. They also show that the peak magnitude appears earlier in time for both media functions f1 and f2. Comparing Fig 3 to Fig 9(B) and Fig 10(A), we see that, while simply having media reports during an epidemic can greatly reduce the severity of an epidemic, it is further mitigated under the optimal reporting intensity. Moreover, in such a scenario, a stronger media report intensity u2* for media function f2 gives rise to a greater reduction in peak magnitude and the total number of infected individuals than the media report intensity u1*.

Discussion

It is known that media reports can play an important role in generating public awareness and promoting disease mitigation measures. Quantifying and evaluating the media impact on the control of emerging infectious diseases is quite challenging. Our study here included the intensity of media reports as a separate compartment, and a modified transmission rate that is reduced using a media function that is bounded between 0 and 1. We also considered a novel media function which depends both on the number of infected individuals I and the intensity of mass media M (where mass media reports were assumed to depend on the reporting of newly infected individuals, ρσE).

We calculated the basic reproduction number of system (1) which is the same as that for the classical SEIR model. This result illustrates that the mass media has no effect on the basic reproduction number, which agrees with the previous studies [11, 24, 27]. We then investigated the threshold dynamics of the proposed model with a general media function effect. We theoretically investigated the optimal control problem by seeking an optimal media reporting intensity to minimize the total infected individuals and the cost of media reporting. The optimal media reporting intensity obtained here indicated that during the early stage of epidemic we should quickly enhance the media reporting intensity, keep it at the maximum level during the time period around the peak of the epidemic, and then decrease the intensity after the epidemic vanishes significantly.

By fitting our proposed model to laboratory-confirmed case data from the 8th Hospital of Xi’an over the first 19 days of the 2009 H1N1 influenza pandemic (to ensure that the data lie in the exponential growth phase of the epidemic), we estimated the unknown model parameters and the basic reproduction numbers without media impact and with two special media functions. We found that the basic reproduction number may be underestimated if the media impact is not considered. We also illustrated that the peak magnitude of the endemic would greatly decrease when the mass media function is considered, which has also been demonstrated in. Sensitivity analysis indicated that the severity of the disease outbreak is sensitive to the parameters associated with the media impact (the reporting rate ρ, media waning rate δ, weight of media effect sensitive to infected individuals α1, and weight of the media effect sensitive to media items α2) besides the epidemiological parameters (like transmission rate β and the recovery rate γ). In particular, the outbreak severity is more sensitive to the weight α1 than weight α2 (shown in Figs 5–7), indicating that more response to the number of infected individuals will lead to greater reductions in the peak magnitude and the total number of infections.

For the particular media functions f1 and f2, we observed no obvious differences in the epidemic curve (shown in Fig 3) when the optimal reporting rate was not considered, though from the point of R-square or AIC for the Least-Square case, the model with f2 can fit the observed data better than the model with f1. However, with the optimal media reporting rate the optimal epidemic curves are quite different. In particular, the shape of the optimal control u1* and u2* are similar, but u2* is greater than u1*, causing the optimal solution of system (4) with f2 to be less than the optimal solution of system (4) with f1. This means that there is optimal media reporting function such that the number of infected individuals and costs reach the minimum, which helps design an optimal news releasing patterns that mostly affect individuals’ behaviour changes, and hence result in the infection significantly decline.

In previous studies, media campaigns have been characterized as a dynamic variable. Media data can be collected to inform the mass media compartment [8, 25]. In our current study, we include a separate compartment for mass media reports. We will consider the incorporation of mass media data in future work.

In summary, we extend the classical SEIR model by incorporating the media as a separate compartment and through the modification of the transmission rate by a media factor associated with not only the number of infected individuals I but also the media items M. Through the inclusion of the media compartment and the modified transmission rate βf(I, M) we can use model (1) to study the effects of I and M separately. In this study, we focus our work on understanding the media impact on the transmission of 2009 H1N1 in Shaanxi, China, and explore the efficiency of optimal control on the media reporting rate. Ultimately, we find that response to the number of infected individuals I will lead to greater reductions in the peak magnitude and the total number of infections. We also find that the optimal media reporting intensity should be enhanced early in the outbreak, be kept at the maximum level during the time period around the peak of the epidemic, and then be decreased after the epidemic reaches a low level to ensure a minimal level of infection in a population.