Epidemics on Networks with Large Initial Conditions or Changing Structure

In this paper we extend previous work deriving dynamic equations governing infectious disease spread on networks. The previous work has implicitly assumed that the disease is initialized by an infinitesimally small proportion of the population. Our modifications allow us to account for an arbitrarily large initial proportion infected. This helps resolve an apparent paradox in earlier work whereby the number of susceptible individuals could increase if too many individuals were initially infected. It also helps explain an apparent small deviation that has been observed between simulation and theory. An advantage of this modification is that it allows us to account for changes in the structure or behavior of the population during the epidemic.


Introduction
The mathematical study of infectious disease spread has contributed significantly to our ability to design effective interventions to reduce disease spread. Most of the earliest models were based on the assumption that disease transmission occurs as a Poisson process and each transmission reaches an individual chosen randomly from the population. This implicitly assumes that partnership duration is very brief. These models have been modified to account for a number of different effects, such as demographic groups [1].
More recently, attempts have been made to incorporate the ''network'' structure of the population (see, e.g., [2]). Typically these focus on trying to understand the role played by ''highdegree'' individuals (those individuals with many contacts), and they come in one of two flavors: they either continue the assumption of fleeting partnerships (the timescale of individual transmissions is long compared to the timescale of individual partnerships) [1,[3][4][5][6], or they take the opposite limit in which the partnership network is static (the timescale of the epidemic is short compared to the timescale of individual partnerships) [7][8][9][10][11][12][13][14]. These two approaches do not address the intermediate regimes.
Recent work has shown that for susceptible-infectious-recovered (SIR) models, it is possible to unify these two approaches with an ''edge-based compartmental model'' (EBCM) that allows partnership duration to range continuously from zero to infinite [15][16][17] [for susceptible-infectious-susceptible (SIS) models, the picture is more complicated, see for example [18]]. The resulting models are low-dimensional and contain many standard models as special cases [17]. Unfortunately, these models are derived under the assumption that the initial proportion infected is infinitesimally small (while the absolute number infected is sufficiently large that the dynamics are deterministic). It is assumed that by the time the equations are used, any early transients have died away. A consequence of this assumption is that the models break down if R 0 v1 (that is, if the average number of infections caused by an infected individual early in the epidemic is less than 1) or if the initial proportion infected is not negligible.
The failure if the initial proportion infected is not negligible was observed by [14]. This paper used an early (static network) version of the equations of [15] from [7] and compared them with simulation. A small discrepancy in final sizes was noted. This discrepancy was not present for equations of [19], a system requiring O(M) equations where M is the maximum degree, or for another system introduced in [14] which required O(M 2 ) equations.
In the remainder of this paper, we generalize the EBCM equations for the spread of infectious disease assuming a finite proportion of the population is initially infected. We test the resulting equations against simulations, analyze their predictions for different disease scenarios, and investigate the cause of the discrepancies found in previous work. For simplicity we focus on the static network limit. The method we introduce is straightforward to adapt to dynamic networks.

Analysis
We modify the approach of [15] which assumed an infinitesimal initial proportion infected. We adapt the approach to consider a wide range of possible initial conditions. We assume that the dynamics of the epidemic may be treated as deterministic, which means we assume the population is very large and the initial number infected is large enough for the epidemic to behave deterministically. The assumption that behavior is deterministic may be understood qualitatively as equivalent to the claim that no individual has a large enough effect to alter the dynamics of the disease at the population-scale: the time (or even if) a single given individual becomes infected has a negligible impact on the proportion infected. If stochastic effects are still important but R 0 w1, then these equations may become accurate at a later time once sufficient numbers are infected.
We assume the population consists of N&1 individuals. Each is assigned a degree k (independently of degrees of other individuals) with probability P(k) where P defines a probability distribution on the non-negative integers. The network is wired together using the ''Configuration Model'' (or ''Molloy-Reed'') approach [11,20]: each individual is assigned a number of stubs (or half-edges) equal to its degree. Pairs of stubs are then wired together to form edges/ partnerships. It is likely that this algorithm produces a handful of self-loops or repeated edges, but although they may be present, their density (i.e., the probability a given individual is involved in a self-loop or repeated edge) goes to zero like 1=N.
We define a test individual u to be a random individual chosen at time 0. Because we assume that the spread is deterministic, this means that the probability u is in a given state is equal to the proportion of the population in that state. So we focus on calculating the probability u is susceptible, infected, or recovered. We modify u so that it does not transmit to any of its partners if ever infected. This assumption does not affect the probability u is in any given state, but it does prevent a correlation between the statuses of different partners which would be caused by infection traveling through u. This allows us to treat the partners of u as independent and so each partner of u may independently transmit to u. It is important to note that this assumption has no impact on the probability u is in any given state and therefore, it does not affect our calculation of the proportion of the population in each state. Further discussion of the test individual is in [21].

Variables and Parameters
We introduce our variables, our parameters, and their definitions in table 2. The starting point is the test individual u. The remaining variables and parameters can be broadly divided into four groups. N S, I, and R denote the proportion of the population in each state, or equivalently the probability that the test individual u is in each state.
N h, w S , w I , and w R give the probability a partner of u has a given status and the probability the partner has transmitted to u: h is the probability the partner has not transmitted to u and w S , w I , and w R give the probability the partner has not transmitted to u and is susceptible, infected, or recovered respectively. N P(k) tells us the probability a random individual has degree k, while S(k,0) tells us the probability a random individual of degree k is initially susceptible. The function y(x)~P S(k,0) P(k)x k encodes P and S. We define SKT~P k kP(k) to be the average degree.
N We have two disease parameters to consider: b, the transmission rate, and c, the recovery rate.
Given our definitions, y(h(t)) is the probability u is susceptible at time t. By noting that y(h(t))~S(t), we will be able to close our system of equations.
The main distinction between this approach and the previous EBCM approach [15] is that we use just the initially susceptible individuals to define y while the earlier work took y(x)~P P(k)x k , the probability generating function for the degree distribution. The earlier work then assumed the disease had already been spreading prior to time 0 and defined h to be the probability that a partner has never transmitted [so h(0)v1] whereas here we take h to be the probability that a partner has not transmitted given that it had not prior to time 0 [so h(0)~1].

Equation Derivation
We will find a closed system of equations based on these variables. We begin by looking at S(t). If the test individual u has degree k and is susceptible at t~0, then the probability it is susceptible at some later time is h(t) k . If we do not know k or whether u is susceptible at t~0, then the probability u is susceptible at time t is the sum over all k of the product of the probability u is initially susceptible S(k,0) with the probability u is still susceptible h k . We have S(t)~P k S(k,0)P(k)h(t) k~y (h(t)).

S(t)~y(h(t))
We know that R solves _ R R~cI. We also have a conservation rule that SzIzR~1, so I~1{S{R. Thus our equations are Assuming h(t) is known, then this system with an initial condition for R completely defines S, I, and R. This is shown in the flow diagram in figure 1. Other formulations are possible, for example S~y(h), _ I I~{y'(h) _ h h{cI, _ R R~cI. However we find our system to be preferable because it minimizes the number of differential equations.
In order to close this system of equations we need an equation giving h. Recall that h(t) is the probability a partner v of u that had not yet transmitted to u by time 0 has still not transmitted by time t. This is broken into three disjoint sub-compartments h~w S zw I zw R based on whether v has not transmitted and is susceptible, infected, or recovered. Because w I is the probability v has not transmitted and is infected, it is straightforward to see that _ h h~{bw I . So if we can find w I in terms of h, then we arrive at a single equation for h, which can be used to provide h for the S, I, and R equations.
To do this, we use the fact that w I~h {w S {w R and find w S and w R in terms of h. We turn to figure 2. The recovery rate is c and the transmission rate is b, so we have _ w w R~{ c _ h h=b. We can integrate this, and using the fact that h(0)~1 we find w R~c (1{h)=bzw R (0). To find w S in terms of h, we note that the probability u has an edge to a node that is susceptible at time t~0 is w S (0). The probability the susceptible partner has degree k is kP(k)S(k,0)= P kP(k)S(k,0), so the probability an initially susceptible partner is susceptible at some later time is (1). We arrive at with h(0)~1. This completes our system. Our final closed system of equations is where h(0)~1, and R(0) is given by the initial conditions. These equations lead to earlier equations of [7,15,22] if R 0 w1 and 1{w S (0), w I (0), w R (0) and 1{h(0) are all infinitesimally small. If R 0 v1, the error caused by the approximations w R (0)~0 and w S (0)~1 is comparable to the actual number infected.
Final size relation. The final size relation assuming small initial condition is well-known [11]. The final size relation for larger initial conditions has recently been found [21] in a more general case not assuming constant transmission and recovery rates. It can be derived easily for this model by setting _ h h~0. We find

Model Validation
We now compare our model with simulations for populations that satisfy the Configuration Model/Molloy-Reed model assumptions. Although an earlier version of these equations was found to have minor discrepancies [14], we show that once we appropriately account for the initial condition, the calculation becomes correct.
Final Size Comparison. To show that our new equations accurately calculate the impact of the initial conditions, we first consider epidemic spread in networks with the same degree distribution as in [14] (table 1), but with varying numbers initially infected and varying population sizes. We then consider the impact of selecting high or low degree nodes as the earliest infected individuals, using networks whose degree distributions more clearly show the impact of biased selection of the initial individuals.
We run a large number of simulations for each number of initial infections. For each simulation we generate a new network. Our simulation technique is similar to those recently described by [8,19,23]. In the Configuration Model framework, each node is assigned a degree, nodes are given stubs (or half-edges), and then stubs are randomly paired together. In the simulations we use, each node is assigned a degree, nodes are given stubs, and then the disease begins to spread in the network before stubs are paired. Each time the disease transmits along a stub that stub is joined to a randomly selected as-yet-unpaired stub. If the partner is susceptible, then it becomes infected. If not, nothing happens. Once stubs are paired they remain in their edge. This approach is equivalent to constructing the network in advance and then following the disease, but it is more efficient computationally because it only constructs those parts of the network the disease traces.
Randomly selected initial infections. We first consider varying numbers of randomly chosen initial infected individuals. In figure 3 we take the degree distribution from table 1.
We randomly select a proportion r of the population to initially infect. We have S(k,0)~1{r for all k, so y(x)~(1{r) P k P(k)x k . Similarly we have w S (0)~1{r. Because the epidemic begins with no recovered individuals, we take w R (0) ls0. We take b~0:1 and c~0:2 (though all that matters for the final size is their ratio).
We take populations of 100, 1000, and 10000 and perform many simulations. We compare the final sizes observed with the final size relation of equations (3) and (4). The equations are derived in the infinite population limit, but in figure 3 we see that even with populations of only 100 they give a good prediction of the observed behavior. As the population size increases, the distribution becomes narrower and the simulations collapse more tightly around the prediction.
Biased initial infections. To show that the approach we have derived can also be applied to cases where the initial infected individuals are selectively chosen based on their degree, we use a different degree distribution that helps highlight the effect. We take P(1)~P(9)~1=2. We consider two options. In the first approach, individuals with higher degree are preferentially selected. To do the selection, we choose an individual with probability proportional to the square if its degree, and infect it. We repeat this until a proportion r of the population is infected. In the second approach individuals are chosen with probability proportional to the square of their inverse degree until a proportion r is infected. We take b~0:1 and c~0:6.
Using these rules, we clearly see that S(k,0) is not uniform. Instead, for the case where individuals are selected with probability proportional to their squared degree, we find that S(k,0)~a k 2 where a solves P k P(k)a k 2~1 {r: We find w S (0)~P kS(k,0)P(k)a k 2 = P k kP(k). In the case where individuals are selected with probability inversely proportional to their squared degree, we find that S(k,0)~a 1=k 2 where a solves P k P(k)a 1=k 2~1 {r, and w S (0)P k S(k,0)P(k)a 1=k 2 = P k k P(k). We compare predictions and simulations in populations of 1000 individuals in figure 4. In the limit of a negligible initial proportion infected, the final size of epidemics in these networks is about 4%. As we increase the number of initially infected individuals, we increase the final size because of these individuals and because of the additional infections they lead to. At small amounts, increasing the initial number of high degree nodes has a much larger impact on the final size because they cause more additional infections. However, as the amount of infection initially present is increased this effect becomes less important: the high degree individuals would become infected anyway. So the largest gain in final size comes from infecting low degree individuals who would not receive an infection from their partners. The ''kinks'' that occur just above 50% initially infected are because effectively all individuals of high (left) or low (right) degree are initially infected.
Dynamic Calculation. We now look at the performance of the dynamic equations. The dynamic prediction is more easily affected by noise than the final size prediction, so we use larger population sizes. We again take the degree distribution of table 1. We begin with 5% infected, either randomly chosen, or chosen as before proportional to the square of the degree. A comparison of individual simulations with theoretical predictions is in figure 5. The theory accurately predicts the dynamics of epidemics in large populations, but there is significant stochasticity in smaller  Table 2. The variables we need to calculate the epidemic dynamics. In all of these u is a test individual: randomly chosen from the population and modified so that it cannot infect others, although it can become infected.

Variable/parameter Definition
Test Individual u A randomly member of the population chosen at time t~0 who is prevented from transmitting to its partners.

S(t)
The proportion of the entire population that is susceptible.

I(t)
The proportion of the entire population that is infected.

R(t)
The proportion of the entire population that is recovered.
The probability a random partner v of u that did not transmit to u by time 0 has not transmitted to u by time t.
The probability a random partner v that did not transmit to u by time 0 is susceptible at time t.
The probability a random partner v that did not transmit to u by time 0 is infected at time t but has not transmitted to u.
The probability a random partner v that did not transmit to u by time 0 is recovered at time t and never transmitted to u.

P(k)
The probability an individual has degree k.

SKT~P k kP(k)
The average degree.

S(k,0)
The probability an individual with degree k is initially susceptible.
y(h(t))~P k S(k,0)P(k)h(t) k The probability that the test individual u is susceptible at time t. In a large population this should equal S(t).

Intervention Impact
We can use our equations to compare the impact of several interventions. We consider an epidemic spreading in the population, and at some intermediate time we introduce a change in the disease or population. Because the system changes at a time with a non-negligible amount of infection in the population, the equations derived assuming a negligible proportion infected fail.
We consider three interventions that may be introduced at time t 1 . All are aimed at ''halving'' the transmission rate, but they do this in different ways. In mass-action based models, these would all have the same effect. We can clearly identify differences using our approach.
N An intervention that reduces b by a factor of 2. N An intervention that reduces b so that per-contact transmission probability b=(bzc) is reduced by a factor of 2.
N An intervention that eliminates half of the partnerships randomly.
The distinction between the first two comes from the fact that partnerships have duration, so reducing b by a factor of 2 does not reduce the infection probability by a factor of 2. The expected number of transmissions an individual sends to a given partner is b=c, but once the partner is infected, the subsequent transmissions are irrelevant. If we use mass action assumptions however, each transmission is to a replacement partner. In each case, halving b halves the total number of transmissions, but when partnerships have nonzero duration, a larger proportion of transmissions have no effect. The probability of transmitting at least once along a static partnership is b=(bzc). So to reduce infection probability by a given factor requires a larger reduction to b. Note that the work of [24] suggests that in Configuration Model networks the final size of our second and third intervention will be the same, (but that in clustered networks it will be different).
We will demonstrate our approach in all three cases, restarting the calculations when the intervention is put into place. In all cases, this allows us to use the conditions at t 1 to predict the final size. We take y 0 (x)~y(x), h 0 , w S,0 , w I,0 , and w R,0 to correspond to time less than t 1 . We use a subscript of 1 for times after t 1 . We Figure 3. Results of simulations for N~100, 1000, and 10000 individuals. The solid curve gives our prediction for the final sizes of epidemic in a large population. Colors are log scale giving probability of that particular epidemic size. Each simulation is for a new network generated using the P(k) from table 1, with b~0:1 and c~0:2. We randomly select a proportion r of the population to initially infect and compare final size with the prediction of theory. The number of simulations for each r for N~100, 1000, and 10000 was 50000, 11500, and 3000 respectively. To show that this is sufficient to resolve the distribution, for :175vrv:225 there were 2000000, 50000 and 10000 simulations performed respectively for each N. This only slightly improves the tails of the distribution. Note that when the initial number (not proportion) of infections is small, a large fraction of simulations end without an epidemic. doi:10.1371/journal.pone.0101421.g003 Case 1. We begin by reducing b by a factor of 2 at time t 1 . Until time t 1 , we are solving the original equations. By solving the original system until t 1 we have h 0 (t 1 ). The probability an individual of degree k is susceptible at time t 1 is S(k,t 1 )~S(k,0)h 0 (t 1 ) k . So our new y(x) is y 1 (x)~P k P(k) S(k,0)h 0 (t 1 ) k x k~y (h 0 (t 1 )x). We take our new h 1 to have h 1 (t 1 )~1. The intervention we are doing has no impact on the probability a partner is in any given state. w S , w I , and w R keep the same proportion, but are scaled up to sum to 1 so each is scaled by h 0 (t 1 ). For example, w S,1 (t 1 )~w S,0 (t 1 )=h 0 (t 1 )~w S,0 (0)y' 0 (h 0 (t 1 ))=y' 0 (1)h 0 (t 1 ).
We restart the solutions with these new values. Case 2. The total probability of transmitting to a partner is b=(bzc). For this intervention we change b so that b=(bzc) is reduced by a factor of 2 at time t 1 . This proceeds exactly as above except that the new value of b must be smaller.
Case 3. When we delete half the edges at random, we do not affect the probability that a random partner is in any given state. So the w variables rescale in the same way as for changing b in the previous cases. However, y undergoes a more significant change. As a starting point, consider P 1 (k 1 ), the probability an individual has degree k 1 after edges are deleted. This depends on P 0 (k 0 ), the probability of having k 0 edges prior to deletion. The relation is The probability the individual has degree k 1 and is susceptible at time t 1 is So if we restart the calculations at t~t 1 we have S(k 1 ,t 1 )~Q(k 1 )=P 1 (k 1 ) and (in general if we delete edges with probability p and keep with probability q~1{p, then the new function is y 1 (x)ỹ 0 (½h 0 (t 1 )½pzqx)). Using this new y 1 (x), the same system of equations holds. Figure 6 compares these strategies. As anticipated, the final sizes resulting from cases 2 and 3 are identical, regardless of the time of intervention (indeed this is straightforward to show from the final size relation). However, we see that the dynamics are significantly different.
There are other ways we could capture these interventions mathematically. We note that in both case 1 and case 2, we could also capture the intervention by simply using a step change in b. Case 3 could be captured by reducing w S , w I and w R each by half and placing the other half into a new inactive compartment w X for the deleted edges. However, each of these is ad hoc and dependent on the precise details of the case. Using the approach presented above, we have a standard approach that will apply across a wide range of interventions.

Bifurcation analysis
We try to gain a better understanding of the epidemic threshold and what happens to the final size as the initial proportion infected is increased. Consider the final size relation found from To study what happens when w I (0)w0 (but possibly arbitrarily small), we begin by first analyzing the structure of the dynamical equation for h under the assumption that w S (0)~1 and w R (0)~0, taking hv1. These assumptions contradict our initial conditions, but understanding the bifurcation in this system first will lead to an easier understanding of the full system with w I (0)w0. If w S (0)~1 and w R (0)~0, the differential equation for h becomes Clearly h~1 is an equilibrium. Close to h~1, we write h~1zm, so y'(h)=y'(1)~1zmy''(1)=y'(1)zm 2 y'''(1)=2y'(1)z O(m 3 ). The value of y''(1)=y'(1) plays a key role in this analysis. Its biological interpretation comes from looking at a random infected individual's partner v early in the epidemic. If v is infected by that random infected individual, then y''(1)=y'(1) is the expected number of additional partners v has (its excess degree).
Substituting into the equation for the equilibrium we have So there is a bifurcation as the bracketed term passes through zero, when (bzc)=b~y''(1)=y' (1). This is the well-known epidemic threshold [11]. The bifurcation is transcritical and corresponds to R 0 increasing through 1. If y''(1)=y'(1)v (bzc)=b (that is R 0 v1) then m is positive and the corresponding equilibrium has hw1 and is unstable, while the equilibrium at h~1 is stable. In our case, we will not observe hw1 because h is a probability. If however y''(1)=y'(1)w(bzc)=b (that is R 0 w1), then the corresponding equilibrium has hv1 and is stable while the equilibrium at h~1 is unstable.
Previous studies found these approximate fixed points for h, and used a slightly different definition for h such that h(0) was slightly less than 1. So the biologically implausible prediction is that for R 0 v1 the value of h increases towards the stable fixed point at 1. For R 0 w1 the prediction is more meaningful: h decreases away from the unstable fixed point at 1 to the lower, stable fixed point. When we do a more careful analysis, we now have h(0)~1 and w I (0)w0. The stable fixed point that was at 1 for R 0 v1 will be slightly decreased to h 0 v1~h(0). So h will decrease to h 0 when R 0 v1. For R 0 w1 the unstable fixed point that was at h 0~1 is slightly increased, so h 0 w1~h(0). So h will decrease away from h 0 . When the number of infections decays immediately (R 0 v1) we care about a small error in the prediction caused by the fact that the initial proportion infected is not exactly zero, while if the number of infections grows (R 0 w1), the error caused by treating the initial proportion infected as asymptotically 0 is insignificant.
We are now able to consider the effect of realistic initial conditions. We keep w R (0)~0, but take w I (0) to be a small positive number with h(0)~1 and w S (0)~1{w I (0). The bifurcation diagram changes slightly. Compared to the equations assuming w S (0)~1, this has the effect of decreasing _ h h slightly, so the equilibrium values shift as shown in figure 7. In fact, strictly speaking, there is no bifurcation for w I (0)w0.
Below the bifurcation value of m, the equilibrium h~1 is slightly reduced to a h 0 v1, but remains stable. The solution with initial condition h~1 converges to this equilibrium. Above the bifurcation, the h~1 equilibrium is slightly increased to a h 0 w1 and is unstable. The other equilibrium with smaller h is stable and its location is decreased slightly. Thus the solution with initial condition h~1 decreases and converges to the stable solution.

Results and Discussion
In this paper we have extended previous mathematical methods for epidemic spread in networks to allow for simple models that can capture the impact of a non-negligible initial condition.
Our system of equations (1) and (2) are mathematically simple and can be solved numerically with standard tools. Changes in the population's degree distribution affect y, but do not otherwise alter the structure of the equations, and in particular, the population can have arbitrarily large maximum degree without requiring any increase in the number of equations.
Our modeling approach accurately predicts the size and dynamics of simulated epidemics with arbitrary sized initial conditions. The approach allows us to compare interventions introduced during the epidemic. We further see that the discrepancy found by [14] apparently results from the fact that the initial proportion infected was nonnegligible.
We are able to investigate the dynamical structure of the equilibria found in the equations, showing how the epidemic threshold is modified when a non-negligible proportion of the population is infected. This helps resolve the apparent contradiction in earlier work in which h, a measure of the probability of having escaped transmission, could increase in time. On biological grounds h must be monotonically decreasing. The mathematical inconsistency resulted from the assumption of a small initial condition. For R 0 v1, in the small initial condition limit, there is an attracting equilibrium value where the cumulative amount of infection is very small. For a reduced initial infection, the equilibrium moves, going to no infection in the limit. The previous analyses effectively froze this equilibrium at its asymptotic limit and then considered a small, but nonzero initial amount of infection. So the initial condition was on the wrong side of the equilibrium and the mathematical model attempts to return to its frozen equilibrium of no cumulative infection rather than approaching the true equilibrium with a small number of recovered individuals. When R 0 w1, this effect does not matter because the frozen equilibrium is repelling and still on the correct side of the initial condition.
One of the most obvious applications of these results is to the understanding of the impact of an intervention that begins after a disease has established itself. This has been a weakness of network models for some time: the earliest low-dimensional models could only calculate static quantities such as the final size of epidemics assuming no intervention, while more recent approaches that calculate the dynamics [7,15,22] have been restricted to the assumption of asymptotically small initial conditions, again with no change in the population behavior. Because we now have a lowdimensional model that can account for large initial conditions, we can use this to restart our calculations when an intervention is to be implemented, or we can use the final size relation to quickly compare intervention effectiveness.
We have analyzed the bifurcation structure of the final size relation, and used this to explain an apparent discrepancy in earlier work if R 0 v1. The previous models that assumed small initial condition also implicitly assume that R 0 w1. This resulted in a disturbing prediction for R 0 v1 that transmissions could be reversed as time progresses, and infected individuals are uninfected. Once we correctly account for the initial condition this apparent discrepancy disappears.
If variance is large enough, then there may be a small number of very high degree individuals who have a macroscopic effect on the dynamics. Usually, increasing the population size will ''drown out'' the effects of individuals. However, when the variance is sufficiently large increasing the population size results in a small number of much higher degree individuals. These again have a macroscopic effect on the dynamics. Deterministic predictions will not be accurate: for example, how long the highest degree individual remains infected will influence the final size. The work of [8] rigorously studied the equations using small initial conditions, and showed that if all moments up to the fifth moment were finite, then these equations are accurate in the limit of a large network. More recently [25] has shown that the equations remain valid with much weaker assumptions. The equations will fail if the degree distribution is too broad exactly because these very rare individuals are able to have a population-scale impact on the dynamics, so the specific time and duration of their infections matter.
We expect to be able to recover from this challenge however. After the epidemic has run for a short period of time in a Configuration Model network, all of these high degree individuals have been infected and recovered. The remaining population will have significantly reduced moments. At this point, the stochastic effects are ''frozen in'': the dynamics are now deterministic. We can use the observed conditions at this time to initialize our new system of equations.
There are some assumptions implicit in our derivation that deserve further attention. The model fails if h(0), w S (0), w I (0), or w R (0) depend on degree of u. So if for example, we select high degree individuals and then infect their partners (leaving the high degree individuals uninfected), the model will not account for the fact that higher degree individuals are more likely to have infected partners at t~0. The approach will fail. This is discussed in more detail in [26].
To be clear, the model does not fail if the initial individuals infected have higher (or lower) degree. This simply affects the initial conditions. Indeed, we expect that if the infection is initially spreading stochastically in the population, and we set t~0 to be when enough cases are infected to have deterministic behavior, we will see that at t~0 a disproportionate number of higher degree individuals have been infected. So w I and w R may initially be larger than I and R.
We finally note that in [15,16], a number of generalizations of the EBCM equations were considered. The basic approach we have used here can be applied to any of those generalizations.

Author Contributions
Conceived and designed the experiments: JCM. Performed the experiments: JCM. Analyzed the data: JCM. Contributed reagents/materials/ analysis tools: JCM. Contributed to the writing of the manuscript: JCM.