Modeling and Controlling Epidemic Outbreaks: The Role of Population Size, Model Heterogeneity and Fast Response in the Case of Measles

Abstract

Modelers typically use detailed simulation models and vary the fraction vaccinated to study outbreak control. However, there is currently no guidance for modelers on how much detail (i.e., heterogeneity) is necessary and how large a population to simulate. We provide theoretical and numerical guidance for those decisions and also analyze the benefit of a faster public health response through a stochastic simulation model in the case of measles in the United States. Theoretically, we prove that the outbreak size converges as the simulation population increases and that the outbreaks are slightly larger with a heterogeneous community structure. We find that the simulated outbreak size is not sensitive to the size of the simulated population beyond a certain size. We also observe that in case of an outbreak, a faster public health response provides benefits similar to increased vaccination. Insights from this study can inform the control and elimination measures of the ongoing coronavirus disease (COVID-19) as measles has shown to have a similar structure to COVID-19.

1. Introduction

The ongoing novel coronavirus disease 2019 (COVID-19) pandemic has shown the importance of efforts to control outbreaks before they become a pandemic. A potential COVID-19 vaccine provides hope for containing the pandemic. Public health officials need to plan for the implementation of vaccination programs while also considering the vaccine hesitancy (necessity, safety and freedom of choice) [1]. Insights from other well-known diseases, especially ones with a similar epidemic progression structure, can help public officials in decision making. Measles is one such example. Due to the structural and progression similarities, researchers have suggested understanding and utilization of origin, models and vaccines of measles for controlling the ongoing COVID-19 pandemic [2,3,4,5,6,7,8,9,10,11,12]. Researchers have shown that there are similarities between the protein S in COVID-19 with protein F in measles [3,4]. Many have suggested using the measles, mumps and rubella (MMR) vaccine as an intervention [5,6,7,8,9,10,11]. They develop epidemic models following the same standard structure as measles where the population is divided into susceptibles, exposed, infectious and recovered groups [12]. Research suggests that the occurrence of COVID-19 might follow a similar pattern to recent measles outbreaks in the United States [13]. For this reason, in this paper we study the model and control measures of epidemic outbreaks in the case of measles to transfer the knowledge to COVID-19.

Vaccination coverage gaps have resulted in spikes in epidemic outbreaks of preventable infectious diseases among which measles has become the most alarming worldwide [14]. Measles is a highly contagious disease, which has a secondary attack rate of more than 90% among susceptible people [15]. It can have severe consequences, especially in young infants and pregnant women, such as encephalitis and death. For every 1000 reported cases of measles, approximately one case of encephalitis [16] and two to three deaths may occur [17]. The introduction of the measles, mumps and rubella (MMR) vaccine decreased the incidence of measles in the 1990 s and measles was declared eliminated in the United States in 2000 [18]. However, since then the number of measles cases in the U.S. has increased to levels last seen in 1992 [19]. Vaccine hesitancy combined with large scale measles outbreaks in Europe have increasingly exposed susceptible American travelers, reintroducing measles to the United States dozens of times each year by infected travelers and leading to a resurgence of measles in the U.S. [19,20,21,22,23]. Despite the continued efforts on awareness and vaccination campaigns, vaccine hesitancy remains a challenge in certain population pockets for religious, medical or personal belief reasons [23,24,25,26]. Fortunately, other key public health response activities such as quarantining of exposed individuals and exclusion of susceptible contacts from specific high-risk settings have shown the potential to help in outbreak control in addition to the vaccination [26,27].

The analyses of outbreak notifications and response data are a key part of an outbreak response [28] and modeling can help decision makers analyze different response strategies. Complex modeling is time intensive and thus often does not support the quick and ad-hoc decision making needed for an outbreak response. This paper analyzes three issues to assist modelers and policy makers with the recent outbreaks, which are two modeling issues (the necessary model size and complexity) and one public health issue (the benefit of a faster response to epidemic outbreaks). To describe the issue of model size, we note that the median outbreak size is five individuals [22] although the largest outbreaks in the U.S. infected a couple of hundred [29,30,31]. This raises the question whether an outbreak simulation needs to include 100, 1000 or 100,000 individuals. For a more detailed agent based simulation, e.g., [32], the related question is what to include in the simulation (i.e., where to draw the boundary of the study population). For example, when simulating an outbreak initiated by a high schooler, the question would be whether to only simulate the individuals in the child’s school or whether to include the elementary school of the initial teen’s younger sister in the simulation. For our stochastic model of epidemic outbreaks, we prove that the outbreak size converges in a distribution as the model population increases and then use the simulation to determine that a population size of 1000 is more than enough for measles.

The most related stream of research on model size is on the critical community size, the minimum size within which measles can persist. Early research suggested that the critical community size was 250,000–500,000 [33,34,35]. Keeling and Grenfell [36] examined outbreaks in 60 towns in England and Wales from 1944 to 1968 to arrive at a critical community size for measles of 250,000–400,000. However, this research refers to the pre-vaccination era where significant fractions of children became infected each year [37] (an estimated number of 3 to 4 million people per year in the United States [38]) instead of the isolated outbreaks we see today. We now turn to the issue of model complexity or heterogeneity. Measles is a well-studied disease and most models incorporate some form of heterogeneity such as age-structure, spatial heterogeneity, vaccination coverage and heterogeneity for the transmission probability between people [32,39,40,41,42,43,44,45]. The cost of including many and detailed forms of heterogeneity are additional to model complexity and simulation run-time. To capture the range of potential heterogeneity we use a stylized model that divides the population into two subgroups and varies the mixing rates between the subgroups and within each subgroup, an approach similar to that taken by one of the earliest measles models with heterogeneity [39] and a later study on the role of heterogeneity [32]. We compare this with a model with a homogeneous, undivided population. We prove that outbreaks are longer when the population is divided into two subgroups and that heterogeneity increases infection risk. However, using simulation, we find that the actual effect of varying the model heterogeneity is quite small.

Turning to the public health interventions for measles, the effectiveness of the measles vaccination in school and university settings [46,47,48,49] and vaccination campaigns [50] is well-studied. The increased rates of religious or philosophical exemptions to vaccination requirements [30,31] make it important to consider alternative interventions. However, studies of other interventions such as better case detection or a faster public health response are missing. Our study includes three intervention variables, which are the vaccinated fraction, the probability of detecting and reporting a case and the number of days’ delay in the public health response. Our numerical simulations show that a faster public health response can be just as important as increasing vaccination coverage.

2. Materials and Methods

We built a model to determine the size of the outbreak, when it was detected and how various interventions or model characteristics affected it. We modeled the progression of measles using an SEIR (Susceptible-Exposed-Infectious-Recovered) model with a 10-day exposed (or incubation) period when susceptible individuals came into contact with the infection but were not able to spread it to others and a subsequent eight day infectious period when they were able to spread the infection to other susceptible individuals [51,52]. We modeled measles transmission using a stochastic simulation model. Every day there was a chance that the outbreak could be detected and contained shortly thereafter with a vaccination program (or closure of the school). At the start of the simulation there was a single infectious individual and we focused on the final size C of the outbreak. We assumed a constant population of size n, of which a fraction v was unvaccinated. We examined the effect of community structure on measles transmission and how people connected and communicated. To study the extremes of community structure we compared the case of homogenous mixing (where everybody was part of the same community group) and the case where the population was divided into two equal-sized community groups. In the latter case, an individual was more likely to infect another individual in the same community group than one in the other group, where this relative likelihood was the mixing parameter, $\mathsf{\phi }$. This was similar to the non-random mixing from Glasser et al. [53] where the inter-school contacts and disease specific infection probabilities decreased based on distance in the case of measles. When $\mathsf{\phi }=0$ there was no mixing between the two community groups and when $\mathsf{\phi }=1$ there was no preferential mixing between the groups (there was effectively only one big community group). However, we kept the basic reproductive number, $R_0$, the same in all cases. We also considered the less stylized case of transmission between two schools, a smaller elementary school of 500 students and a high school with 2000 students.

We assumed each infectious individual had a daily probability q of being detected and reported. Throughout the paper we use detecting and reporting a case interchangeably. Note that sometimes the outbreak was never detected due to underreporting [54] and not detecting the imported cases [55]. Under homogenous mixing (a single community group), the simulation was halted $ϵ$ days after the first case was detected and reported. With two community groups, we assumed that public health officials separately stopped the spread in each group. Specifically, the spread in each group was halted $ϵ$ days after the first detected case in that group and $\sigma$ = 5 days after the first detected case in the other group, whichever occurred first.

2.1. Model

Here we describe the equations that describe our model. Below is also a table of notation. We start by describing the equations for the case with a homogenous population. We let n be the size of our population and v the fraction unvaccinated. Thus, before the outbreak there were vn susceptible individuals in the population. On day t, this population was divided into S(t) susceptible individuals, E_k(t) individuals in day $k\in \left\{1,2,\dots ,10\right\}$ of the incubation period, I_k(t) individuals in day $k\in \left\{1,2,\dots ,8\right\}$ of the infectious period and R(t) individuals who had recovered from an infection. Note that $S\left(t\right)+{{\displaystyle \sum }}_{k=1}^{10}{E}_k\left(t\right)+{{\displaystyle \sum }}_{k=1}^8{I}_k\left(t\right)+R\left(t\right)=vn$. We denoted the total number of infectious individuals on day t as T(t). On day 1, we had exactly one infected individual who was in the first day of the incubation period. The system dynamics are then given by Equations (1)–(9), where N(t) was the number of new infections on day t.

$$T\left(t\right)={\displaystyle \sum }_{k=1}^8{I}_k\left(t\right).$$

(1)

$$S\left(t\right)+{\displaystyle \sum }_{k=1}^{10}{E}_k\left(t\right)+T\left(t\right)+R\left(t\right)=vn.$$

(2)

$$S\left(1\right)=nv-1,E_1\left(1\right)=1.$$

(3)

$$S\left(t\right)=S\left(t-1\right)-N\left(t\right).$$

(4)

$$E_1\left(t\right)=N\left(t\right).$$

(5)

$$E_j\left(t\right)=E_{j-1}\left(t-1\right),\mathrm{for}j\in \left\{2,3\dots ,10\right\}.$$

(6)

$$I_1\left(t\right)=E_{10}\left(t-1\right).$$

(7)

$$I_k\left(t\right)=I_{k-1}\left(t-1\right),\mathrm{for}k\in \left\{2,3,\dots ,8\right\}.$$

(8)

$$R\left(t\right)=R\left(t-1\right)+I_8\left(t-1\right).$$

(9)

Equations (10) and (11) describe the number of new infections N(t) on day t. These were distributed binomially with each susceptible individual on day t − 1 having a probability $p\left(t\right)$ of being infected. We calculated $p\left(t\right)$ as the total number of infectious individuals, $T\left(t-1\right)$ times the basic reproductive number R₀ and divided by the length of the infectious period, 8, and the size of the population. This way, the first infectious individual created R₀ secondary cases in expectation if the entire population was susceptible.

$$N\left(t\right)~\mathrm{Binomial}\left(S\left(t-1\right),p\left(t\right)\right).$$

(10)

$$p\left(t\right)=\frac{T\left(t-1\right)\left(\frac{R_0}{8}\right)}{n}.$$

(11)

Equations (12)–(14) describe the final size of the outbreak, C. Here, Z(t) was the number of cases reported on day t. We then found the time, t, which gave us the first reported case of the outbreak ($\min \left\{t:Z\left(t\right)\ge 1\right\}).$ Thus, on day τ, ε days after the first reported case, public health authorities stopped the spread of the outbreak by coordinating various activities such as press releases about the index case, isolation, school closures, polymerase chain reaction (PCR) tests and postexposure prophylaxis (PEP) administration [56]. The coordination of these activities generally takes time [56], thus we used a lag term, ε, to account for the delay. The lesser the value of ε meant the faster public health response. The outbreak size C was then the number of infected and removed individuals on that day. The number of reported cases per day, $Z\left(t\right)$, was distributed binomially with each infectious individual having a probability q of being detected and reported each day.

$$Z\left(t\right)~\mathrm{Binomial}\left(T\left(t\right),\mathrm{q}\right).$$

(12)

$$\tau =ϵ+\min \left\{t:Z\left(t\right)\ge 1\right\}.$$

(13)

$$C=nv-S\left(\tau \right).$$

(14)

Now we describe the small modifications to the model we needed to make for the case where the population was divided into two subgroups. These are summarized in Equations (15)–(22). We let $n^x,{\mathrm{S}}^{\mathrm{x}},E_k^x,I_k^x,R^x,T^x,N^x,Z^x,{\tau }^x,C^x$ and $p^x$ be the respective variables for subgroup x = 1 or x = 2. Initially, the only infected individual was in subgroup 1 and was in the first day of the incubation period. Recall that the mixing parameter φ was the ratio of the probability of infecting a susceptible in the other subgroup to the probability of infecting a susceptible in the same subgroup. Thus, we adjusted the infection probabilities p^x(t) by considering both infections from the same subgroup and the other subgroup. In addition, we not only had a public health intervention stopping the spread in a subgroup ε days after the first reported case in that subgroup but also $\sigma$ days after the first reported case in the other subgroup.

Table 1. Notation.

Population, n

Fraction unvaccinated, v

Susceptibles, S(t)

Individuals in day k of the incubation period, Ek(t)

Individuals in day k of the infectious period, Ik(t)

Recovered individuals, R(t)

Total number of infectious individuals, T(t)

New infections, N(t)

Probability a susceptible individual is infected, p(t)

Basic reproduction number, $R_0$

Number of reported individuals, Z(t)

Probability an infectious case is reported each day, q

Lag between the first reported case in same subgroup and stopping, ε

Stopping time, τ

Mixing parameter, $\phi$

Lag between the first reported case in the other subgroup and stopping, $\sigma$

below summarizes the notation used in this manuscript.

$$S^1\left(1\right)=n^{1}v-1,E_1^1\left(1\right)=1.$$

(15)

$$S^2\left(1\right)=n^{2}v.$$

(16)

$$n=n^1+n^2.$$

(17)

$$p^1\left(t\right)=T^1\left(t-1\right)\left(\frac{R_0}{8}\right)\frac{1}{n_1}\frac{1}{1+\phi }+T^2\left(t-1\right)\left(\frac{R_0}{8}\right)\frac{1}{n_1}\frac{\phi }{1+\phi }.$$

(18)

$$p^2\left(t\right)=T^2\left(t-1\right)\left(\frac{R_0}{8}\right)\frac{1}{n_2}\frac{1}{1+\phi }+T^1\left(t-1\right)\left(\frac{R_0}{8}\right)\frac{1}{n_2}\frac{\phi }{1+\phi }.$$

(19)

$${\tau }^1=\min \{ϵ+\min \left\{t:Z^1\left(t\right)\ge 1\right\},\sigma +\min \left\{t:Z^2\left(t\right)\ge 1\right\}\}.$$

(20)

$${\tau }^2=\min \{ϵ+\min \left\{t:Z^2\left(t\right)\ge 1\right\},\sigma +\min \left\{t:Z^1\left(t\right)\ge 1\right\}\}.$$

(21)

$$C=C^1+C^2.$$

(22)

2.2. Simulations

We conducted numerical experiments to verify our analytical results and provide additional insights. For our simulations, we considered a stylized case with a focus mostly in the U.S. We utilized the average elementary school size for the subgroup size to obtain the total population size [57]. By using the average vaccination level of 90%, we obtained the fraction unvaccinated [58]. Despite being a mandatory reported disease, public health authorities still struggle vastly with underreporting. We gathered daily reporting probabilities based on studies about the underreporting of measles in other developed countries [54,59,60]. We found the delay in the reporting and the public health action from two case studies of measles in Colorado, U.S. [56]. For the baseline of the mixing parameter we chose 5, the middle value of the range, which was [0, 1] and for the lag between the first reported case in the other subgroup and stopping, we chose five days.

Table 2. Parameters.

Parameter	Value	Source
Population, n (people)	1000	[57]
Basic reproduction number, $R_0$	18	[37]
Fraction unvaccinated, v (%)	10%	[58]
Incubation period (days)	10	[51]
Infectious period (days)	8	[52]
Mixing parameter, $\phi$ (%)	0.5
Probability of the case being reported each contagious day, q (%)	0.1	[54,59,60]
Lag between the first reported case in the subgroup stopping, ε (days)	3	[56]
Lag between the first reported case in the other subgroup and stopping, σ (days)	5

presents the summary of the parameter values.

We tracked the number of infected cases, C, and the number of infected cases in the second community group, the one that was initially not infected, C₂, when we assumed two community groups. Specifically, we examined the expected size of the outbreak, E[C], the probability the outbreak was at least five, P [C ≥ 5] and the probability that the outbreak spread beyond the initial community group, P [C₂ ≥ 1]. In the post elimination era, the median number of outbreak cases is five; thus, five was the threshold outbreak size we used in our metric [22].

We also examined the impact of public health interventions such as vaccinated fraction of the population, daily reporting probability and faster public health response. We varied the vaccination coverage using the different vaccination rates across states [61]. The values of daily reporting probability was changed based on the range in different underreporting studies [54,59,60]. In order to measure the impact of faster response, we investigated results by reducing the delay term from 3 days to 0 days. During the sensitivity analysis, all parameters not varied in a figure or table were at their baseline value as given in Table 2.

We simulated several scenarios to conduct a sensitivity analysis for the less stylized case of transmission between two schools, a smaller elementary school of 500 students and a high school with 2000 students. We assumed that initially there was a single infected case in the elementary school as measles cases were more likely to appear in the elementary school [62].

3. Results

3.1. Theoretical Results

In this section, we analyze the impact of population size on the model and the model heterogeneity. For the population size, we show that the epidemic outcome measures converged to finite limits as the population size increased. For the model heterogeneity, we show that the outbreaks were larger and longer in the inhomogeneous case.

Recall that $T=T\left(t,n\right)$ denoted the infectious population and $N\left(t,n\right)$ the number of new infections, where we may drop the dependence on t and n for convenience. We let $C\left(t,n\right)$ denote the size of the outbreak at time t, the number who were infected or recovered at this time. We let $Y\left(i\right)$ denote whether individual i became infected at time t so that the number of new infections at time t was $N={{\displaystyle \sum }}_{i=1}^{nv-C}Y\left(i\right)$. For given T, the $Y\left(i\right)$ was independent and identically distributed with $Y\left(1\right)~\mathrm{Bernoulli}\left(T\beta /n\right)$, where $\beta$ was the effective contact rate.

3.1.1. Convergence as Population Size Increased

Our goal in this subsection is to prove for the homogenous case that the infection process converged to a finite limit as the population size became large. For notational convenience, let $C^{\prime }=C\left(t+1,n\right)$ and note that $C\left(1,n\right)=1$ and $C^{\prime }=C+N$.

Theorem 1.

For all t, $T\left(t,n\right)$$N\left(t,n\right)$and$C\left(t,n\right)$converge in distribution as$n\to \infty .$

Lemma 1.

For all$t,$$E\left[C\left(t,n\right)\right]=O\left(1\right)$as$n\to \infty .$If$E\left[C\right]=O\left(1\right)$, then$E\left[C^{\prime }\right]=O\left(1\right)$as$n\to \infty .$

Proof. 

Since $C\left(1,n\right)=1,$ it suffices to show that if $E\left[C\right]=O\left(1\right)$ then $E\left[C^{\prime }\right]=O\left(1\right)$ as $n\to \infty$ and induct on t. Since $N\le {{\displaystyle \sum }}_{i=1}^{nv}Y\left(i\right)$ as $E\left[N|T\right]\le nvE\left[Y\left(1\right)|T\right]=\beta vT$ a.s. Hence, $E\left[N\right]=E\left[E\left[N|T\right]\right]\le \beta vE\left[T\right].$ Therefore, $E\left[C^{\prime }\right]=E\left[C\right]+E\left[N\right]\le E\left[C\right]+\beta vE\left[T\right].$ Since $T\le C$ a.s., $E\left[T\right]\le E\left[C\right]$ and $E\left[C^{\prime }\right]\le E\left[C\right]\left(1+\beta v\right)$, thus proving the induction hypothesis. □

Lemma 2.

Suppose that$X\left(n\right)\ge 0$a.s. and$E\left[X\left(n\right)\right]\to 0.$Then$X\left(n\right)\to 0$ in distribution.

Proof. 

Since $X\left(n\right)\ge 0$ as $X\left(n\right)=\left|X\left(n\right)-0\right|$ a.s. Hence, $E\left[X\left(n\right)\right]=E\left[\left|X\left(n\right)-0\right|\right]\to 0.$ Finally, convergence in mean implies convergence in distribution. □

Lemma 3.

If$w\left(n\right)\to w\left(\infty \right)$and$y\left(n\right)\to 0$, then${\left(1+w\left(n\right)/n\right)}^{n\left(v-y\left(n\right)\right)}=e^{vw\left(\infty \right)}.$

Proof. 

Note that $\log \left({\left(1+w\left(n\right)/n\right)}^{n\left(v-y\left(n\right)\right)}\right)=\left(v-y\left(n\right)\right)n\log \left(1+w\left(n\right)/n\right)$ equals $\left(v-y\left(n\right)\right)n\frac{w\left(n\right)}{n}\left(1+o\left(1\right)\right)=\left(v-y\left(n\right)\right)w\left(n\right)\left(1+o\left(1\right)\right).$ Since this is a product of limits, it converges to $vw\left(\infty \right)$. The claim follows from the continuity of exp(). □

Lemma 4.

Suppose$T\left(t,n\right)\to T\left(t,\infty \right)$in distribution. Then$N\left(t,n\right)\to \mathrm{Poisson}\left(T\left(t,\infty \right)\beta v\right)$ in distribution.

Proof. 

The characteristic function of a random variable $X$ is $E\left[e^{isX}\right]$ as a function of $s.$ Let ${\phi }_n\left(s\right)=E\left[e^{isN\left(t,n\right)}\right]$ be the characteristic function of $N\left(t,n\right).$ Define the characteristic function of $\mathrm{Poisson}\left(T\left(t,\infty \right)\beta v\right)$ as

$$\begin{gathered} \phi \left(s\right)=E\left[e^{is\mathrm{Poisson}\left(T\left(t,\infty \right)\beta v\right)}\right]=E[E[e^{is\mathrm{Poisson}\left(T\left(t,\infty \right)\beta v\right)}|T\left(t,\infty \right)]] \\ =E\left[e^{\left(T\left(t,\infty \right)\beta v\left(e^{is}-1\right)\right)}\right] \end{gathered}$$

where we use the fact that the characteristic function of $\mathrm{Poisson}\left(\lambda \right)$ is $e^{\left(\lambda \left(e^{is}-1\right)\right)}$. By Levy’s continuity theorem, it suffices to show for all $s>0$, that ${\phi }_n\left(s\right)\to \phi \left(s\right)$ as $n\to \infty$.

Note that $e^{isN}={{\displaystyle \prod }}_{i=1}^{nv-C}{e}^{isY\left(i\right)}$ a.s. Thus,

$$\begin{gathered} E\left[e^{isN}|T,C\right]=E{\left[e^{isY\left(1\right)}|T,C\right]}^{nv-C}={\left(1-\frac{T\beta }{n}+\frac{T\beta }{n}{e}^{is}\right)}^{n\left(v-C/n\right)} \\ ={\left(1+\frac{T\beta \left(e^{is}-1\right)}{n}\right)}^{n\left(v-C/n\right)} \end{gathered}$$

where we use the fact that the characteristic function of Bernoulli(p) is $(1-p+p{e}^{is}$). By Lemmas 1 and 2, $C/n\to 0$ in distribution. We now construct a joint probability space where $T\left(t,n\right)\to T\left(t,\infty \right)$ a.s. and $C\left(t,n\right)/n\to 0$ a.s. Then, by Lemma 3, $E\left[e^{isN}|T,C\right]\to e^{\left(T\left(t,\infty \right)\beta v\left(e^{is}-1\right)\right)}$ a.s. Since $e^{isN}$ is bounded, we can apply the bounded convergence theorem to take the expectation of both sides and prove the claim. □

Proof of Theorem 1. 

Since $C\left(1,n\right)=1$ and $T$ and $C$ are the sum of previous new infections, we can apply induction using Lemma 4. □

3.1.2. Outbreaks Are Longer with Two Subgroups

We now turn our attention model heterogeneity. There were two reasons why heterogeneity increased the expected outbreak size. The first reason was when a community was divided into two subgroups: the outbreak lasted longer when the community was divided into two subgroups. This was because if a case was detected in one subgroup, public health officials did not intervene in both subgroups at the same time. Instead, they intervened later in the subgroup that did not detect the first case. In the homogenous case, the intervention when the first case was detected affected more people immediately (i.e., the entire population) than in the two subgroup case. To formalize this, let ${\tau }^{\left(1\right)}=\max \left\{{\tau }^1,{\tau }^2\right\}$ and ${\tau }^{\left(2\right)}=\min \left\{{\tau }^1,{\tau }^2\right\}$. Also, define the case detection times in the homogenous case and the case with two subgroups, $\tilde{\tau }=\min \left\{t:Z\left(t\right)\ge 1\right\}$ and ${\tilde{\tau }}^x=\min \left\{t:Z^x\left(t\right)\ge 1\right\},$ respectively.

Theorem 2.

Suppose that delay between detection and intervention is larger if the case is detected in the other subgroup,$\sigma \ge ϵ$and that everything else, the total number of infections and detected cases is the same as in the homogenous case,$Z^1\left(t\right)+Z^2\left(t\right)=Z\left(t\right)$for all t. The stopping time in the homogenous case is then the same as in one of the subgroups but before the stopping time in the other subgroup,$\tau ={\tau }^{\left(2\right)}\le {\tau }^{\left(1\right)}\le \tau +\sigma -ϵ$.

Proof. 

Note that $\tilde{\tau }=\min \left\{{\tilde{\tau }}^1,{\tilde{\tau }}^2\right\}$. Hence, $\tau =\tilde{\tau }+ϵ=\min \left\{{\tilde{\tau }}_1,{\tilde{\tau }}_2\right\}+ϵ={\tau }^{\left(2\right)}$. Obviously, ${\tau }^{\left(2\right)}\le {\tau }^{\left(1\right)}$. The last inequality, ${\tau }^{\left(1\right)}\le \tau +\sigma -ϵ$ follows from the definition of the stopping time (20) and (21). □

Additional intuition for $\tilde{\tau }\le {\tilde{\tau }}^1,{\tilde{\tau }}^2$ (a consequence of $\tilde{\tau }=\min \left\{{\tilde{\tau }}^1,{\tilde{\tau }}^2\right\}$) is the following. The number of person-days of infected individuals required until a case was detected is a geometric random variable with mean 1/q. In populations with more infected individuals (which larger populations tend to have), one had a higher probability of detecting a case and detection occurred more quickly.

3.1.3. Heterogeneity in Infection Risk

We now focus on the second way that model heterogeneity led to larger outbreaks: the variation in the infection risk between individuals in the case with two subgroups. We abstracted beyond two subgroups and considered the general case where the infection probability varied among individuals, $Y\left(i\right)~\mathrm{Bernoulli}\left(p_i\right).$ To ensure a valid comparison, we used the average infection probability in the homogenous case, $p^{*}=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{p}_i=T\beta /n$, ensuring that $E\left[N\right]={{\displaystyle \sum }}_{i=1}^{nv-C}{p}_i=T\beta \left(v-C/n\right)$ in any case. Formally, we again assumed that the $Y\left(i\right)$ were independent given T. Our argument that outbreaks grew more slowly in the homogenous case has two parts (we show later that outbreaks were also more likely to fizzle out in the homogenous case). First, we will show that the variance in the number of new infections, $\mathrm{Var}\left[N\right]$, was maximized in the homogenous case. Second, we go on to show that increases in the variance of the number of infectious individuals, $\mathrm{Var}\left[T\right]$, decreased the expected number of new infections. While we prove both parts below, we were unable to combine them to argue about the outbreak size at any particular time, $C\left(t\right)$, due to difficulty with the induction step.

Lemma 5.

The homogenous case maximizes$\mathrm{Var}[N|T,C]$.

Proof. 

Let $f\left(x\right)=x\left(1-x\right).$ Note that $\mathrm{Var}\left[N|T,C\right]=\mathrm{Var}\left[{{\displaystyle \sum }}_{i=1}^{nv-C}Y\left(i\right)|T,C\right]={{\displaystyle \sum }}_{i=1}^{nv-C}f\left(p_i\right)$. Since $f$ is concave, $\left(1/k\right){{\displaystyle \sum }}_{i=1}^{k}f\left(p_i\right)\le f\left(p^{*}\right)$ where $p^{*}$ is the average infection probability as defined above. Thus, ${\displaystyle \sum }_{i=1}^{nv-C}f\left(p_i\right)\le \left(nv-C\right)f\left(p^{*}\right)=\mathrm{Var}[\mathrm{Binomial}\left(n,p^{*}\right)|T,C]$, proving the claim. □

Lemma 6.

Let$C=T+R$. Then, $E[N|R]$decreases in$\mathrm{Var}\left[T|R\right]$.

Proof. 

By definition,

$$E\left[N|R\right]=E\left[{{\displaystyle \sum }}_{i=1}^{nv-T-R}Y\left(i\right)|R\right]=E\left[E\left[{{\displaystyle \sum }}_{i=1}^{nv-T-R}Y\left(i\right)|T\right]|R\right]=E\left[{{\displaystyle \sum }}_{i=1}^{nv-T-R}{p}_i|R\right].$$

Since $Y\left(i\right)$ are iid, $\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{p}_i=T\beta /n$, $E\left[N|R\right]=E\left[\left(nv-T-R\right)\left(\frac{T\beta }{n}\right)|R\right]$ and thus

$$\begin{array}{ll}E[N|R]\hfill & =E\left[T|R\right]\beta \left(v-\frac{R}{n}\right)-E\left[T^2|R\right]\frac{\beta }{n}\hfill \\ & =E\left[T|R\right]\beta \left(v-\frac{R}{n}\right)-\left(E{\left[T|R\right]}^2+\mathrm{Var}\left[T|R\right]\right)\frac{\beta }{n}\end{array}$$

proving the claim. □

In addition to the expected number of new infections being smaller in the homogenous case, we also showed that the probability of no new infections was larger in that case. This is particularly important at the beginning of an outbreak, which can fizzle out if the first few infected individuals do not manage to infect any others.

Lemma 7.

The homogenous case maximizes$P[N=0|T,C]$.

Proof. 

The proof is similar to that of Lemma 5. Note that $P\left[N=0|T,C\right]={{\displaystyle \prod }}_{i=1}^{nv-C}P[Y\left(i\right)=0|T]={{\displaystyle \prod }}_{i=1}^{nv-C}\left(1-p_i\right)$. Since $\log \left(1-x\right)$ is a concave function, $\left(1/k\right){{\displaystyle \sum }}_{i=1}^k\log \left(1-p_{\mathrm{i}}\right)\le \log \left(1-p^{*}\right)$. Hence, $\log P\left[N=0|T,C\right]\le \left(nv-C\right)\log \left(1-p^{*}\right)=\log P\left(\mathrm{Binomial}\left(n,p^{*}\right)=0\right)$, proving the claim. □

3.2. Numerical Results

3.2.1. The Impact of Population Size

In this section, we not only confirm the convergence proof in Section 3.1 but also provide how large a population size is large enough to be close to converging. Figure 1 examines the effect of changing the population size on the average outbreak size and the probability of a large outbreak, one where five or more become infected. After increasing initially, both outcome measures varied little as the population increased beyond a thousand. This implies that there was no gain to simulating measles outbreaks with larger populations.

3.2.2. The Impact of Heterogeneity

In this section, we confirm that heterogeneity resulted in more infections and numerical results allowed us to determine whether the increase was significant. Figure 1 also shows how the average outbreak size and the probability of a large outbreak was greater for a community divided into two subgroups than a homogenous community while keeping the total population size fixed. Figure 2 focuses on this aspect of the community structure by comparing the distribution of the outbreak size in the two cases. It also shows that large outbreaks were slightly more frequent in the inhomogeneous case. There were multiple reasons including the difference in stopping times between subgroups and the heterogeneity of the force of infection.

Figure 3 further explores the community structure in the case of the two subgroups by varying the mixing parameter, which specified how separate the two subgroups were. The mixing parameter was the ratio of the likelihood that an individual infected someone in the other subgroup rather than someone in the same subgroup. When the mixing parameter was zero there was no mixing between subgroups. In that case, we focused on the subgroup with the initial infection (because it would never spread to the other subgroup) and viewed it as a single homogenous population of half the original size. This partly explains why the expected number of infected cases was lower with no mixing. When the mixing parameter was one, there was complete mixing. This was like a single homogenous population except that, as discussed in the previous figure, there was an additional lag until the outbreak was contained. Figure 3 shows that as the mixing parameter increased, the expected number of outbreak size seemed to slightly increase.

3.2.3. The Impact of Public Health Interventions

We now turn our sensitivity analysis to important parameters for public health interventions. Figure 4 focuses on the effect of changing the vaccinated fraction of the population. For each percentage point that the vaccinated fraction increased, the expected size of the outbreak decreased by approximately 0.35 and the probability of a large outbreak decreased by approximately 3.3 percentage points.

Two other important parameters for public health interventions were the daily probability of detecting an infected case and the lag between the first detected case and a public health response stopping the spread of the outbreak.

Table 3. Additional sensitivity analysis. Here C₂ is the number of infected cases in the other subgroup (not the one with the initial infection). The number of replications was 1000. The standard error was less than 0.18 for E[C] and 0.02 for P [C ≥ 5] and P [C₂ ≥ 1].

		Two Community Groups			One Community
		E[C]	$\mathit{P}\left[\mathit{C}\ge 5\right]$	$\mathit{P}\left[{\mathit{C}}_2\ge 1\right]$	E[C]	$\mathit{P}\left[\mathit{C}\ge 5\right]$
Daily reporting probability, q	0.05	5.9	0.465	0.650	5.2	0.403
	0.10	3.9	0.287	0.554	3.4	0.257
	0.15	3.2	0.199	0.533	2.8	0.159
	0.20	2.7	0.116	0.499	2.6	0.108
Lag between first case reported and simulation stopped,ε	0	3.1	0.212	0.558	2.8	0.185
	1	3.4	0.247	0.531	3.0	0.192
	2	3.7	0.262	0.567	3.4	0.239
	3	3.7	0.251	0.551	3.5	0.242

examines the sensitivity to these parameters for the average outbreak size (E[C]), the probability of a large outbreak with five or more infected (P[C ≥ 5]) and, for the case of a population with two subgroups, the probability the infection spreads from one subgroup to the other (P[C₂ ≥ 1]). While the benefits of increasing the detection probability were significant, we saw diminishing returns for greater increases. Decreasing the lag between detection and response also decreased the outcome measures (except for the probability of the second subgroup being infected, which changed little) but only by a small amount. By the time a case was detected, the outbreak had already spread to the other subgroup. Thus, a few days of delay did not make much difference and the change in P[C₂ ≥ 1] was even less.

We simulated several scenarios to conduct a sensitivity analysis for the less stylized case of transmission between two schools, a smaller elementary school of 500 students and a high school with 2000 students. We assumed that initially there was a single infected case in the elementary school. In

Table 4. Sensitivity analysis for two unequal-sized groups, an elementary school of 500 and a high school of 2000 individuals. The initial infection was in the elementary school and C₂ is the number of cases in the high school. There were 1000 replications. The standard error was less than 0.23 for E[C] and 0.02 for P[C ≥ 5] and P[C₂ ≥ 1].

		E[C]	$\mathit{P}\left(\mathit{C}\ge 5\right)$	$\mathit{P}\left({\mathit{C}}_2\ge 1\right)$
Fraction vaccinated, %, (1-v)	80	8.3	0.621	0.823
	85	5.9	0.464	0.717
	90	4.0	0.308	0.579
	95	2.1	0.076	0.306
Mixing parameter,$\mathbf{\phi }$	0	3.4	0.237	0.000
	0.25	3.8	0.285	0.451
	0.5	3.9	0.297	0.570
	0.75	3.9	0.279	0.659
	1	4.1	0.315	0.677
Daily reporting probability, q	0.05	5.9	0.460	0.639
	0.1	3.8	0.276	0.552
	0.15	3.1	0.200	0.525
	0.2	2.8	0.128	0.492
Lag between first case reported and simulation stopped,ε	0	3.1	0.211	0.537
	1	3.6	0.261	0.572
	2	3.8	0.286	0.556
	3	3.9	0.307	0.572

, we examine different choices for the fraction vaccinated, the mixing parameter, the daily probability of detecting a case and the lag between the first reported case and stopping the outbreak. The outcome measures were the same as in Table 1, with P[C₂ ≥ 1] now denoting the probability that the outbreak spread from the elementary school to the high school. The results had the same trend as those in the case of a population with two equal-sized subgroups we discussed previously. As the vaccinated fraction of the population, the mixing parameter or the response lag increased, the expected number of infected cases, the probability of having at least five infected cases and the probability of having at least one infected case in the high school when the public health authorities took action decreased.

4. Discussion

We presented a stochastic simulation model of a measles outbreak to investigate the impact of the simulated population size, model complexity and public health interventions. We studied this for both homogeneous mixing (i.e., a single community) and an inhomogeneous case where the population was divided into two community groups. We observed the final number of infected cases, the number of infected cases in the second community group where the initial number of infected people was zero and the probability of the infection spreading from the initially infected community group to the second community group. Our theoretical analysis showed that results converged to a finite limit as the population size increased and gave multiple reasons (from the way that outbreaks were stopped to the role of heterogeneity) why the expected outbreak size for a community divided into two subgroups was greater than in the case of a single homogenous community. We also performed a sensitivity analysis on various model parameters as well as suggesting public health interventions such as increasing the probability of reporting a case and decreasing the lag between detecting a case and the public health authority stopping transmission.

We found that the size of the modeled population ceased to matter after exceeding 1000 individuals (Figure 1), matching the theoretical result. This made sense because with 10% unvaccinated, the unvaccinated population was 100 individuals, which was large enough so that the stochastic effects of small numbers were small. For smaller populations, the outbreaks were smaller as they were limited by the population size. Overall, the two subgroup model had slightly larger outbreaks compared with the homogenous model (Figure 2), which also matched the theoretical results. Increasing the mixing parameter in the two subgroup model increased the average outbreak size but the effects were small (Figure 3). In addition, the probability of an outbreak spreading from one subgroup to another was around 50% and was somewhat insensitive to the parameters.

Increasing the vaccinated fraction had a larger effect than increasing the probability of detecting a case, which in turn had a larger effect than decreasing the lag between detection and a public health response in terms of decreasing the size of the outbreak (Figure 4 and Table 3). The asymmetric case with a bigger high school and a smaller elementary school had similar results. Larger outbreaks (i.e., lower vaccination rates or smaller detection probabilities) were also associated with higher probabilities of the infection spreading from one subgroup to another (Table 4).

One limitation of our study was that we only examined two community structures (one with homogenous mixing and one with two subgroups). The hope is that these represent the extremes and other community structures will have outcomes between the two extremes. Our model did not capture the possibility that the unvaccinated individuals were all clustered around the index case in the social network; our model assumed that the unvaccinated were spread randomly within each subgroup. Social network data could evaluate the likelihood of this possibility. This could also be captured in our model by restricting the study population to the smaller set of social contacts of the index case (e.g., maybe 20 individuals with 40% vaccinated). We are not aware of such social network data. Another limitation was that there are no studies to justify a value for the case detection probability in our model. Such a study would be valuable especially because the true extent of measles outbreaks is uncertain (i.e., underreporting) [54] and because we found that increasing the detection probability would be an effective way of decreasing the size of measles outbreaks. In addition, we also did not consider details such as variability in the measles progression across different individuals. However, this was not the focus of our study.

The model and the parameters were mainly based on measles in the U.S. The characteristics of the model especially the values of the parameters can change based on applications to other diseases and geographies. For example, the average MMR vaccination coverage rate is around 90% in the U.S., thus the median number of outbreak size is very low (five people). It might be hard to attain these vaccination coverage numbers in other countries and for other diseases. As a result, the outbreak sizes might be larger.

There are three conclusions we can draw. First, unless the population of interest is actually small, modelers can save a lot of computational effort by limiting their simulated population to 1000 individuals with little cost in accuracy. Second, this is irrespective of the heterogeneity in the population. As heterogeneity increased, the expected results increased slightly. Thus, unless focusing on the impact of specific interactions, one can exclude complex heterogeneity items from models. This would also speed up validation, calibration and sensitivity analysis. Third, in contexts where it is difficult to raise the vaccination coverage (which are many), alternative interventions such as better detection of infections and a faster response to any detected cases should be considered.

COVID-19 has demonstrated the importance of transforming scientific knowledge on epidemic modeling into an effective response that saves as many lives as possible and reduces the financial and health consequences to society and the economy. Many researchers suggest the use of learnings from measles to contain the pandemic as viruses causing both diseases share same protein and have a similar epidemic progression. Hopefully, conclusions drawn from this study can inform policy makers as they make long term plans for the global COVID-19 pandemic. Policy makers and researchers have to consider adjustments. For example, measles is one of the most contagious diseases (with $R_0$ of 18). With a low median size of outbreak, the time it takes for the confirmed cases to be brought to public authorities (in other words, simulation length) is short. For diseases with lower $R_0$ such as COIVD-19, the simulations might need to take longer. Currently, a vaccination for COVID-19 does not exist so the results from this study can be more useful once the vaccination uptake is in further stages such as elimination. This study can also help when considering and evaluating alternative strategies to vaccination during the vaccination uptake. Additionally, insights from this study can help decision makers to deal with future outbreaks efficiently and effectively before they become pandemics.