The Effects of Media Reports on Disease Spread and Important Public Health Measurements

Controlling the spread of influenza to reduce the effects of infection on a population is an important mandate of public health. Mass media reports on an epidemic or pandemic can provide important information to the public, and in turn, can induce positive healthy behaviour practices (i.e., handwashing, social distancing) in the individuals, that will reduce the probability of contracting the disease. Mass media fatigue, however, can dampen these effects. Mathematical models can be used to study the effects of mass media reports on epidemic/pandemic outcomes. In this study we employ a stochastic agent based model to provide a quantification of mass media reports on the variability in important public health measurements. We also include mass media report data compiled by the Global Public Health Intelligence Network, to study the effects of mass media reports in the 2009 H1N1 pandemic. We find that the report rate and the rate at which individuals relax their healthy behaviours (media fatigue) greatly affect the variability in important public health measurements. When the mass media reporting data is included in the model, two peaks of infection result.


Introduction
Controlling the spread of influenza to reduce the effects of infection on a population is an important mandate of public health. Influenza control can be achieved through various strategies including vaccination, the use of drug therapy, hand washing and social distancing-removing oneself as much as possible from the population [1,2]. Vaccines and drug therapies, however, may be ineffective (e.g. resistance to antivirals) or unavailable (i.e. pending development of a new vaccine). Thus, hand washing and social distancing practices (i.e., ranging from moderate to extreme total isolation), which can be used at all times during an infectious disease outbreak are integral in decreasing the probability of contracting and transmitting infection [1].
Mass media campaigns can be used to provide information on current and effective vaccination, drug therapy and social distancing measures [1]. Public health education campaigns, that include informative literature (i.e., pamphlets), posters, newspaper articles and advertisements, radio and television messages, and social media outlets (i.e., twitter, facebook) are used daily to inform the public on current health issues. Mass media outlets can aid in dissemination of this information. Studies of mass media campaigns and healthy behaviour have reported that mass media campaigns can elicit positive behaviour change and even prevent negative behaviour change in individuals [3][4][5][6][7][8][9]. It is therefore concluded that mass media campaigns should be used to inform the public so that behaviour change can result [3][4][5][6][7][8][9][10]. It has been shown that information conveyed by the media is becoming the critical factor as to whether or not a vaccination campaign will succeed [11]. It is believed that the media coverage of the recent SARS and 2009 H1N1 epidemics had an effect on the total spread of these pathogens [4,7,12,13]. More recently, MERS-CoV, Ebola, and H7N9 have been the subject of many media reports.
Pervasive media coverage of social problems may lead to desensitization to media reports [14][15][16]: a diminished emotional responsiveness to a negative or an aversive stimulus after repeated exposure. This same phenomenon can occur with health events -individuals may be more likely to take precautions against becoming ill with the first reports of an epidemic or health threat, but individual sensitivity to disease reports may diminsh with time, and social distancing practices may be relaxed [15,[17][18][19], affecting disease transmission. To better control infectious diseases a better understanding of the effects of mass media on the uptake and waning of social distancing practices is needed.
Mathematical models can be employed to study the effects of media in an epidemic [19][20][21][22][23][24][25][26][27][28][29][30][31][32]. For example, the effects of media reports (i.e., on disease incidence, hospitalizations, deaths) can be incorporated as a function f(S, E, I, R) that directly affects the transmission rate in a susceptible-exposed-infected-recovered (SEIR) epidemiological model [19-27, 32, 33]: This model has provided the basis of the vast majority of past mathematical modelling studies of media during an epidemic. In a previous study on the effects of mass media reports on important public health measures (i.e., peak number of infections, time of peak infection, end of epidemic and total number of infections), we demonstrated that the variability in these measures can vary drastically depending on the media function f(S, E, I, R) chosen [33]. In this study, we then demonstrated that the recommended public health measures to mitigate epidemic outcomes were not consistent over different choices of f(S, E, I, R), and we suggested that models that include a mass media compartment, that better represents population interactions with mass media, should be considered [33]. For example, explicity includes a mass media compartment, allowing a sensitivity analysis over biological assumptions that determine the function f(S, E, I, R). Here, g(M) is a function representing the number of mass media reports M, g(M)S models the interaction that individuals have with media that can change their behaviour, S M is a class of susceptibles that are media aware and are less susceptible to infection (or entirely protected through healthy and social distancing behaviours), and μ represents the effects of media fatigue. A further benefit of Eq (2) is that it enables the inclusion of mass media report data to replace f(S, E, I, R). Studies of the effects of mass media on epidemic outcomes employing systems of equations similar to Eq (2) exist in the literature [28][29][30][31], however, none of these studies have incorporated mass media data into their analysis. Also, quantification of the variability surrounding important public health measures has been ignored.
In this study we extend Eq (2) to include more levels of healthy behaviour, and vaccination. We then study the model in two different ways. First, we employ a stochastic agent-based simulation to study variability in key public health measurements. Second, we incorporate 2009 H1N1 mass media report data from the Global Public Health Intelligence Network (GPHIN) [34,35] into the model and study the effects of this data and media fatigue on the pandemic dynamics in a closed population. From herewithin, we describe all healthy behaviours including isolation, limiting contacts and hand washing as 'social distancing'. Social distancing levels included in the model are chosen to reflect different degrees of protection that healthy behaviour change can provide. Specifically, we include the extremes-no protection and total protection/isolation, and a moderate protection level. The uptake of social distancing behaviours, and vaccination are both assumed to be affected by mass media reports [3,5,6,11].
In the following sections we motivate the model and determine the expected variability in key public health measurements (i.e., peak number of infections, time of peak infection, end of epidemic and total number of infections). A sensitivity analysis of the model is also completed so that parameters that most affect the key public health measurements can be determined. Interestingly, the model parameters representing mass media reports and media fatigue are identified. When mass media data [34] is incorporated into the deterministic model, a twopeak pandemic curve results. This curve has similar qualitative characteristics to that observed over the recent 2009 H1N1 pandemic.

Model
The model considers fully susceptible (S), social distancing susceptible (S 1 ), isolated susceptible (S 2 ), vaccinated (V), exposed (E), infectious (I) and recovered (R) populations. The media (M) is also incorporated as a separate compartment and affects the movement of individuals between the susceptible and vaccinated classes. The model is as follows: A flow diagram of Eq (3) is shown in Fig 1. Briefly, fully susceptible individuals S can be infected (βI) and can take up precautionary measures motivated by interactions (listening, reading, watching) with media reports (M), such as vaccination (ν) or social distancing (α).
Susceptible individuals in the social distancing class S 1 can also be infected, but at a lower rate (β 1 ), and can also take on precautionary measures such as vaccination (ν) and isolation (α 1 ), motivated by exposure to media reports. It is assumed that isolated individuals (S 2 ) are fully isolated and cannot be infected, but can be vaccinated (ν) as a result of exposure to media reports (M). It is assumed that vaccination provides more protection than any level of social distancing, thus, all vaccinated individuals move to the vaccinated class (V). Exposed individuals are produced via infection of susceptible individuals S and S 1 and can become infectious (σ). Infectious individuals can recover from infection (γ). It is assumed that individuals are born (λ) fully susceptible and that social behaviour can wane (q i , i = 1, 2) which moves individuals in higher levels of social distancing to lower ones. It is also assumed that individuals can only die from natural causes (d) (a disease induced death rate is ignored). An increase in media reports can be modelled using various different assumptions. Here, it is assumed that media reports are proportional to a fraction of incoming infectious cases ρσE that are symptomatic and/or reported to a doctor [36]. It is also assumed that the effects of mass media on social distancing behaviours can wane (ρ 1 ) [14,15,37,38]. The relaxation of social distancing behaviours q 1 and q 2 also reflect this behaviour.
Note that Eq (3) ignores covariances and higher order moments i.e. it assumes that all model classes are independent of all others. We include a derivation of Eq (3) in the Appendix that includes higher order moments. For the current study and parameter values, Eq (3) and the model provided in the Appendix do not differ significantly (not shown). Thus, we chose only to present results of Eq (3) here within.
Note that more levels of social distancing can be included in the susceptible populations. Social distancing can also be included in the exposed and infectious classes. For the purposes of the current study we chose to include only three levels of social distancing in the susceptible population-none, moderate, and total isolation. This selection has provided a more analytically tractable model to study mass media effects on social distancing and vaccination behaviours. Furthermore, an optimal number of social distancing levels that would best represent all social distancing behaviours in susceptible, exposed and infected populations is currently unknown. An extension of this model to study a 'continuous social distance' in all model populations, instead of discrete stages like that used here, is a course for future work.

Stochastic model
Eq (3) is deterministic and cannot be used to study variability in infection. We have developed an Agent-based Monte Carlo (ABMC) simulation to study variability in key epidemic measurements when vaccination, social distancing, media and waning behaviour are included. The ABMC simulation models the population at the level of individual agents where each agent contains characteristics associated with infection status, birth, death and social distancing measures. The ABMC simulation moves forward in time using event times: the next time that an individual changes state within the system. To compare to a system of ordinary differential equations (ODEs) exponential distributions for all lifetimes are assumed. Table 1 lists the means of the lifetime distributions, corresponding to the parameters of Eq (3).  Table 1. Briefly, agents in each of the susceptible, exposed, infectious and recovered compartments are assigned event times corresponding to each event that allows an individual from that compartment to change state. To move froward in time the simulation chooses the minimum event time in the population and performs the corresponding event. It then searches for the next minimum event time and so on. The model events and mean event times are listed in Table 1.

Parameter Values
Parameter values for Eq (3) can be found in Table 1. These values agree with the influenza literature and are further informed by the literature surrounding social behaviour responses to media reports.
It is assumed that individuals experience on average 6 days of infection (1/σ + 1/γ) and can transmit the infection for approximately 4 days at the end of the infection time period (1/γ = 4 days). The transmission rate of infection β is calculated using the basic reproductive ratio R 0 where R 0 takes the value 1.5 which lies in the reported range 1.3 to 1.7 [39]. When vaccination is included in the model, it is employed only during the second wave of infection as vaccination was not available during the first wave of the 2009 H1N1 pandemic. It is assumed that on average 25% to 40% of the population is vaccinated [40]. The social distancing uptake and relaxation rate values (α i and q i ) are assumed so that the pandemic curve has a similar magnitude and length to that observed for 2009 H1N1 in Canada [1,2,[39][40][41].
Fig 2 shows data extracted from a study on the 'newsworthiness' of influenza H1N1 as extracted from newspaper homepages across 12 websites from April 29 to May 28, 2009 [16]. The solid line on Fig 2 shows an exponential curve fit to the data points. It is assumed that the media waning rate ρ 1 is equal to the decay rate of this exponential curve fit.
The Government of Canada collects data on emerging public health threats through the Global Public Health Intelligence Network (GPHIN) [34].

Disease free equilibrium (DFE)
The disease free equilibrium (DFE), or uninfected equilibrium, is the equilibrium point at which no infection exists in the system. The DFE for Eq (3) is  Basic Reproductive Ratio The basic reproductive ratio, R 0 , is defined as the number of new infectious individuals produced by one infectious individual when introduced into a fully susceptible population. Generally, if R 0 < 1, a disease will die out, but if R 0 > 1, the pathogen will spread in the population.
There are various methods that can be used to calculate R 0 [42]. Using the survivor function and the next generation methods [42] to calculate R 0 for Eq (3), we obtain where S 0 = λ/d from E 0 . This is an intuitive result as it applies directly to the underlying processes of transmission and progression through the different stages of infection. Here, the first term corresponds to the production of an exposed individual through transmission of infection from an infected individual to a susceptible, and the second term gives the probability that a  [16]. An exponential curve is fit to the data to determine a media fatigue waning rate (red). The resulting equation is f(t) = 0.09911 exp − 0.0469t. Right panel: Media reports. Second wave of an epidemic produced from the fitted value for ρ 1 from [16] and parameter values in [41]. The x-axis is time in days. doi:10.1371/journal.pone.0141423.g002 The Effects of Media on Disease Spread and Public Health Measurements new exposed individual progresses to the infectious compartment. We note here, that the disease free equilibrium E 0 is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1 (see Appendix for proof).

Infected equilibrium
Eq (3) has at least one infected equilibrium if the following equations are satisfied. The Effects of Media on Disease Spread and Public Health Measurements and I and S 1 satisfy Translation to short term epidemics For short term epidemics such as influenza, births (λ) and natural death (d) can be ignored since the total number of births and deaths over the course of the epidemic is very small. In the following sections it is assumed that λ = 0 and d = 0. The DFE then becomes where N is the total population size. There is no infected equilibrium for this reduced system.

Numerical simulations
We considered various scenarios with and without vaccination, social distancing and waning behaviour. Table 2 lists the values of the key epidemic measurements-peak time, peak magnitude, epidemic end, total infected, total vaccinated-resulting from Eq (3) for various combinations of vaccine uptake, media social distancing levels and media waning: (a) social distancing and vaccination are not considered; (b) extending to include one level of social distancing (i.e., α > 0), with no waning in behavioural practices (i.e., q 1 = 0 and ρ 1 = 0) and vaccination (ν = 0); (c) extending to include vaccination (ν > 0); (d) adding an isolated class S 2 (α 1 > 0); (e) adding media media waning (ρ 1 > 0) and movement between social distancing classes through the relaxation of social distancing practices (q, q 1 , q 2 > 0). In general, these studies demonstrated that, as expected, vaccination uptake and social distancing practices will reduce the impact of an epidemic, and that waning behaviour will lessen these benefits.

Variability in infection
Variability in the key epidemic measurements, as calculated from 100 runs of the ABMC simulation are also shown in Table 2. Overall, the results suggest that public health programs should be implemented to consider time intervals of 10-50 days for peak infections, and 20-160 days for the end of the epidemic, depending on the control and social distancing being practiced.
We also see that the variability in the peak magnitude, total number of infections, and total number vaccinated is affected by these practices-decreases and increases in variability (calculated using the coefficient of variation = (standard deviation)/mean) can be observed. It is interesting to note that the addition of public health control measures can increase the variability in the total number of infections and vaccinations. This may be due to the fact that individuals can move between social distancing classes many times during the epidemic, providing greater uncertainty in these outcomes of the epidemic. This will affect public health preparedness strategies. It is important to understand what drives this variability, so that public health programs can focus on certain attributes of an epidemic, so as to minimize the variation.

Sensitivity Analysis
To determine which parameters most affect model the model outcomes, a sensitivity analysis employing Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCC) is performed. These methods are said to be the most efficient methods for performing sensitivity analysis with variation in multiple parameters simultaneously [43]. The outcomes that are of interest for Eq (3) include the peak time of the epidemic (or time to the peak), the end time of the epidemic (or length of the epidemic), the peak magnitude of infectious and exposed individuals, and the total number of infectious individuals. Fig 4 shows the PRCC values for Eq (3) assuming a constant population size N. This figure shows that R 0 is the most important parameter for each outcome. However, the total number of infections over an epidemic is also sensitive to changes in the rates at which individuals become infectious (σ) and recover (γ). The figure also shows that the peak magnitude and the total number of infections are both sensitive to the media parameter ρ, and the epidemic end time is sensitive to the media fatigue rate ρ 1 . Overall, the sensitivity analysis results show that to minimize the effect of the epidemic, public health programs should strive to (1) reduce transmission (effectively reducing R 0 ), (2) effectively use drug therapies that can shorten the infectious period (effectively decreasing σ and increasing γ), (3) increase the number of reported cases ρ, and (4) provide educational programs throughout an epidemic that will help minimize the waning of positive social distancing practices ρ 1 . These results demonstrate that role of mass media can significantly affect the severity and length of an epidemic/pandemic. Thus, attention should be given to studying how best to use mass media during an epidemic, so that epidemic measurements can be minimized.

Mass Media Data and Media Waning
In the previous sections is was assumed that the mass media level of reporting was proportional to the number of report symptomatic cases, related to a fraction of incoming infectious cases ρσE. The structure of Eq (3), however, allows for the incorporation of real mass media data into the model like that compiled by GPHIN (Fig 3) [34]. Fig 5 shows the results of Eq (3) when social distancing and waning behaviour are included in the model with (right column) and without (left column) vaccination and the media data [34]. Here, the English language media report levels are used (see Fig 3, dashed line) (top row), and the media data is then extended to consider various scenarios of media reporting until the end of the pandemic: a constant level equal to the last collected data point (second row), a decay after the final data point (third row), setting the report level to zero at the end of the data set (fourth row), and setting the number of reports to zero at the declared end of the pandemic (week 47 [44]). Note that when vaccination is included in the model (right column), the vaccination rate for each of the scenarios is between 25% to 40% of the population [2], and vaccination is included only in the second wave of infection. Fig 5 demonstrates that for all cases, two waves of the pandemic occur, with the magnitude of the peaks varying by the model assumptions. Interestingly, these dynamics agree with what was observed in many countries globally during the 2009 H1N1 pandemic-two waves of infection with two distinct peaks were observed in this timeframe. For a qualitative comparison we plot the total number of confirmed cases in Canada during the 2009 H1N1 pandemic [44] that shows this two wave behaviour (see Fig 6).  In addition to English, French is another official language in Canada. When the French language reports are included in the model similar results are found-two peaks of infection occur, however, there is a smaller number of total infections due to the fact that there are more media reports used in the media equation and the total population N being considered was unchanged. A further simulation considering the same population size and all media reports in all languages (Fig 3 solid line) results in a a single small wave of infection (not shown). Again, this is an intuitive result since the model assumes that if there are more media reports, more individuals will uptake social distancing practices and vaccination. However, if we increase the total population size, which better represents a population that is affected by media reports in all of the languages tracked by GPHIN, two peaks of infection will result.

Discussion
Mathematical modelling studies of the effects of mass media on an epidemic either include a specific function that decreases the transmission rate or a mass media compartment [19][20][21][22][23][24][25][26][27][28][29][30][31][32]. In a recent study, we demonstrated that epidemic outcomes can vary greatly depending on the mass media function chosen, making it difficult to determine any quantification of the variability in key epidemic measurements for public health [33]. It was concluded that models that include a mass media compartment should thus be employed [33]. Mathematical models that include a mass media compartment also enable the study of media fatigue, and can be further modified to include actual mass media report data. In this study we developed a model that includes a media compartment.
Eq (3) includes a mass media compartment, vaccination and two susceptible classes that practice social distancing at different levels. The model also includes the effects of media fatigue The Effects of Media on Disease Spread and Public Health Measurements and waning social distancing practices. Table 2 demonstrates that social distancing and vaccination benefit a population (a decreased level of infection results). This table also shows, however, that media fatigue and waning behaviour can lessen the positive benefits of these intervention measures (last column).
An agent-based Monte Carlo simulation based on Eq (3) was employed to study the variability in key epidemic measurements that must be considered in public health policies and planning. Results showed that the expected outcome of the simulations is quantitatively similar to the solution of deterministic system of ordinary differential equations of Eq (3), but that, variability in the key epidemic measurements can be quite large (Table 2). This in turn will affect public health resources and planning.
A sensitivity analysis of Eq (3) demonstrated that the parameters associated with mass media reports and media fatigue can significantly affect epidemic outcomes. Thus, more time and effort should be given to the study of the effects of these reports during epidemics/pandemics so that these effects can be accounted for in public health policies and preparedness plans.
In the "Mass Media Data and Media Waning" section we incorporated mass media data, compiled by the Global Public Health Intelligence Network [16,34], into Eq (3). This study demonstrated that two waves of infection can occur when the GPHIN mass media data is included in the model. This is an interesting result for two reasons. First, two waves of infection were experienced by many countries globally during the 2009 H1N1 pandemic, including Canada. Second, previous studies of the two wave behaviour have mainly focused on the effects of school holidays during the pandemic, which effectively reduced transmission and then increased it again a couple months later, producing two waves of infection [44][45][46][47]. Here, we have shown that media fatigue can also produce two waves of infection with similar qualitative dynamics to that observed over the 2009 H1N1 pandemic.
As stated above, models incorporating school closure over the summer months can generate two waves of infection [47]. Differences in mixing between age groups, demography and mobility [48], and now, mass media reports and waning behaviour can also contribute to the two wave phenomenon. It is of interest to determine what the magnitude of these contributions to the double peak behaviour can be. This is a course for future work. This would also lead to a study of the optimal use of media reporting in public health [32,49].
Eq (3) includes two levels of social distancing in the susceptible classes. However, social distancing behaviours can be interpreted to occur over a continuous scale from moderate to intensive levels of social distancing. Also, exposed and infectious individuals can also practice social distancing measures. For example, an infected individual may elect to stay home from work, and distance themselves from family members in the same household as well. To model continuous social distancing in all classes we propose a new model employing partial differential equations that keeps track of all classes of infection over time and social distancing space. Analysis of this model is current work. In future work, the model will be augmented to include asymptomatic infections, where individuals are infectious but do not demonstrate symptoms. This will affect transmission, as social distancing will not be enhanced if infection is asymptomatic, and thus, mass media messages that encourage social distancing will be ever more important.
Mass media messages can relay information that will induce positive or negative social distancing practices. Mass media reports can also affect vaccination uptake. For example, the reporting of the benefits and risks may induce vaccination uptake in some individuals, but induce vaccine hesitancy in others. A study of the effects of the type of mass media message on social distancing behaviour is greatly needed, and is a course for future work.
Finally, in this study we focused our work on understanding the effects of mass media data and media fatigue in the context of the 2009 H1N1 pandemic. Similar studies of other infectious diseases would also be of interest. MERS-CoV and Ebola in the current day are a specific interest.

Stochastic differential equations
As an example, we now derive a system of stochastic differential equations for the standard SEIR model as derived in [50]. We construct a stochastic differential equation model from the continuous time Markov chain (CTMC) model for an SEIR model by considering time as continuous on tε[0, 1), the collection of discrete random variables S(t), E(t), I(t), and R(t) for susceptible, exposed, infectious and recovered classes or states, respectively. Transition probability from state to state is called the infinitesimal transition probability, the change in time is sufficiently small that only one event occurs in the time period, and depends only on the state of the system at the previous time: Markov property. The probability of a transition can be described by Eq (13) From the equations for the infinitesimal transition probabilities, Eq (14), we can derive the forward Kolmogorov differential equations. From the forward Kolmogorov differential equation we derive a probability generating function (pgf), a moment generating function (mgf) and finally the differential equations for the mean and higher order moments. The derivation for the forward Kolmogorov differential equation from Eq (14), representing the rate of change of the transitions from state to state follows. The first step is to derive a discrete equation from Eq (14) for sufficiently small Δt for p i,j,k,l (t + Δt), the transition probability, p i;j;k;l ðt þ DtÞ ¼ lp iÀ1;j;k;l ðtÞDt þ bði þ 1Þkp iþ1;jÀ1;k;l ðtÞDt þdði þ 1Þp iþ1;j;k;l Dt þ sðj þ 1Þp i;jþ1;kÀ1;l ðtÞDt þdðj þ 1Þp i;jþ1;k;l ðtÞDt þ gðk þ 1Þp i;j;kþ1;lÀ1 ðtÞDt þdðk þ 1Þp i;j;kþ1;l ðtÞDt þ dðl þ 1Þp i;j;k;lþ1 ðtÞDt þð1 À ðl þ bik þ di þ sj þ dj þ gk þ dk þ dlÞÞp i;j;k;l ðtÞDt; ð15Þ where p i,j,k,l (t) represents the probability of a susceptible, exposed, infectious and recovered population of size i, j, k, and l, respectively. Next, p i,j,k,l (t) is subtracted from both sides, the expressions are divided by Δt, resulting in Eq (16), p i;j;k;l ðt þ DtÞ À p i;j;k;l ðtÞ Dt ¼ P ½lp iÀ1;j;k;l ðtÞS i E j I k R l þ bði þ 1Þkp iþ1;jÀ1;k;l ðtÞS i E j I k R l þdði þ 1Þp iþ1;j;k;l S i E j I k R l þsðj þ 1Þp i;jþ1;kÀ1;l ðtÞS i E j I k R l þ dðj þ 1Þp i;jþ1;k;l ðtÞS i E j I k R l þgðk þ 1Þp i;j;kþ1;lÀ1 ðtÞS i E j I k R l þ dðk þ 1Þp i;j;kþ1;l ðtÞS i E j I k R l þdðl þ 1Þp i;j;k;lþ1 ðtÞS i E j I k R l Àlp i;j;k;l ðtÞS i E j I k R l À bikp i;j;k;l ðtÞS i E j I k R l Àdip i;j;k;l ðtÞS i E j I k R l À sjp i;j;k;l ðtÞS i E j I k R l Àdjp i;j;k;l ðtÞS i E j I k R l À gkp i;j;k;l ðtÞS i E j I k R l Àdkp i;j;k;l ðtÞS i E j I k R l À dlp i;j;k;l ðtÞS i E j I k R l þ oðDt Dt : Next take sum over (i, j, k, l) and take lim Δt ! 0 and arrive at the forward Kolmogorov differential equation, Eq (17), dp i;j;k;l ðtÞ dt ¼ lp iÀ1;j;k;l ðtÞS i E j I k R l þ bði þ 1Þkp iþ1;jÀ1;k;l ðtÞS i E j I k R l þdði þ 1Þp iþ1;j;k;l S i E j I k R l þsðj þ 1Þp i;jþ1;kÀ1;l ðtÞS i E j I k R l þ dðj þ 1Þp i;jþ1;k;l ðtÞS i E j I k R l þgðk þ 1Þp i;j;kþ1;lÀ1 ðtÞS i E j I k R l þ dðk þ 1Þp i;j;kþ1;l ðtÞS i E j I k R l þdðl þ 1Þp i;j;k;lþ1 ðtÞS i E j I k R l Àlp i;j;k;l ðtÞS i E j I k R l À bikp i;j;k;l ðtÞS i E j I k R l Àdip i;j;k;l ðtÞS i E j I k R l À sjp i;j;k;l ðtÞS i E j I k R l Àdjp i;j;k;l ðtÞS i E j I k R l À gkp i;j;k;l ðtÞS i E j I k R l Àdkp i;j;k;l ðtÞS i E j I k R l À dlp i;j;k;l ðtÞS i E j I k R l : Next, replace the deterministic variables i, j, k, l with the state variables of S, E, I, R to arrive at the probability generating function (pgf) Eq (18) dP dt ¼ lðS À 1ÞPðS; E; I; R; tÞ þ bIðE À SÞ @ 2 P @S@I þ sðI À EÞ @P @E þdð1 À SÞ @P @S þ dð1 À EÞ @P @E þ gðR À IÞ @P @I þdð1 À IÞ @P @I þ dð1 À RÞ @P @R : Here, convert the pgf to the mgf, since the method of derivation of the moment differential equations is simpler from the moment generating function, defined as M(θ, ϕ, ψ, z, t) = P(S, E, I, R, t)e Sθ + Eϕ + Iψ + Rz . From the probability generating function (pgf), take the total derivative of M(θ, ϕ, ψ, z, t) to get function Eq (19).
@M @t ¼ lðe y À 1ÞM þ bðe Ày À 1Þ @ 2 M @y@c þ dðe Ày À 1Þ @M @y þsðe cÀ À 1Þ @M @ þ dðe À À 1Þ @M @ þ gðe zÀc À 1Þ @M @c þdðe Àc À 1Þ @M @c þ dðe Àz À 1Þ @M @z : ð19Þ To derive the system of stochastic differential equations for S, E, I, and R, take the derivative of Eq (19) with respect to the variable that corresponds to the moment required for the differential equation; for example, if E[S] is required, take the partial derivative of Eq (19) w.r.t θ and evaluate the results at (θ, ϕ, ψ, z) = (0,0,0,0), If cov(S, I) = 0, the ordinary differential equation is returned and the state variables of the system are independent. Obviously this cannot be solved explicitly, as it is not closed, differential equations for the higher order moments must be considered, in an attempt to solve the equation. Take the derivative of Eq (19) w.r.t S and I, and end up with It is obvious that this equation contains the third order moment. The system of stochastic differential equations is dE½SI dt ¼ lE½I À bE½IIS À 2dE½SI þ sE½SE À gE½SI dE½EI dt ¼ bE½IIS À sE½E À sE½EI þ sE½EE À 2dE½EI À gE½EI dE½SE dt ¼ lE½E À bE½SI À bE½SEI þ bE½SSI À 2dE½SE À sE½SE: To choose which higher order moment equations to explicitly define, look at which 2 nd order moments appear in the equations of the first order moments. Next look at which equations of 2nd order moments appear in the already existing 2nd order moment equations.