{'fraction_non_alphanumeric': 0.026222932065520187, 'fraction_numerical': 0.006314303755153586, 'mean_word_length': 4.8619638580448505, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling Insights from COVID-19 Incidence Data: Part II -- Why are compartment models so accurate?

The SIR-compartment model is among the simplest models that describe the spread of a disease through a population. The model makes the unrealistic assumption that the population through which the disease is spreading is well-mixed. Although real populations have heterogeneities in contacts not represented in the SIR model, it nevertheless well fits real U.S. state data at multiple points throughout the pandemic. Here we demonstrate mathematically how closely the simple continuous SIR model approximates a model which includes heterogeneous contacts, and provide insight onto how one can interpret parameters gleaned from regression in the context of heterogeneous dynamics.


I. INTRODUCTION
Differences in rates of contact and of susceptibility between individuals strongly affect both their likelihood of catching COVID-19 and their health outcomes once infected. Incorporating these heterogeneities into models of disease spread is essential to understand the differential impacts of the disease upon different subpopulations, such as Black or Brown Americans [1,2], nursing home residents [3], incarcerated [4] and unvaccinated individuals [5]. As well as ensuring a more equitable understanding of the disease, baseline questions about whether, e.g. it was safe to reopen schools in the middle of the pandemic, can not be definitively answered without considering the different levels of vulnerability of the communities affected: students, their caretakers, teachers and school staff.
However, heterogeneity-capturing models contain many unknown parameters that are difficult to fit to real data, and are hard to interpret once fit. Accordingly, public health departments continue to make predictions about the progress of the epidemic and about the effectiveness of social distancing based on so-called well-mixed models (reviewed in [6]), and these models can be made to fit existing case data very well. They suffer, however, from having assumptions that are too simplistic to honestly reflect complex social behaviors inherent in disease spread, and as such, parameter fits should be interpreted with caution. In the well-known SIR model [7], which is the focus of this paper, a group of susceptible individuals (S ) transition via contact with infectious individuals to the infectious (I ) group, and after some time recover or are removed (R). The compartments evolve according to: while dR dt = γI ensures that, neglecting disease and other mortality, the total number of individuals S + I + R = N tot remains constant. The susceptibility coefficient β represents the number of infections caused by a single infected individual in an otherwise susceptible population in unit time. γ is a basic recovery rate. For the above equations to be valid the population should be well-mixed: everyone in the population interacts with everyone else at all times. Given geographic considerations and changes in mixing behavior during a pandemic, this assumption cannot be true. Accordingly, much modern epidemiological modeling, including of COVID-19 [6,8,9], has focused on the role of heterogeneous contacts, either by ramifying compartments or by using networks to model connections between individuals [10][11][12][13][14][15][16][17][18]. Although inclusion of heterogeneities can drastically affect e.g. thresholds for herd immunity [11], the multiplication of parameters that occurs when heterogeneities are added to models makes it hard to validate these predictions against real data. Previous work has shown that under certain conditions heterogeneous models can be approximated by well-mixed models [19][20][21], but there is limited data showing these conditions are met by real epidemics.
In Part I of this paper, we showed that the COVIDcase curves from different US states and during different surges can be clustered into between 4 and 6 groups. The collapse of case curves from different states to a small number of master curves suggests, that in spite of heterogeneous COVID transmission rates and impacts, relatively low model complexity is needed to reproduce the overall growth and the decay in number of infectious cases in a surge. Alignment of case data from different states involved translating data in time, and rescaling number of cases by a population size that was detected during data clustering (Part I), and that we interpreted to be the size of the subpopulation through which COVID-transmission was occurring. It follows that any mathematical model to describe this data collapse, must have the same symmetries; and these symmetries are present within the SIR model [6].
We therefore postulated that the SIR model, described above, might be able to describe some of the families of case dynamics. We fit the SIR model to data from the first COVID surge, specifically the cumulative number of infections detected in each US states from the beginning of their respective outbreaks in February or March of 2020, up to May 20, 2020 (by which time stay at home orders had been relaxed, changing the transmission rate of the disease [22]). In total, the SIR model fits the US state data for 2 of the 4 initial COVID phase clusters, accounting for 26 of the 52 states and territories in the data set (see Fig 1 and Fig 2)   California's cumulative COVID case data (blue) is reasonably approximated by the SIR model (orange) both in its initial stage and during the Omicron surge (left and right panels, respectively). Initial COVID data is highly dependent on testing capacity which was highly variable during the initial wave of COVID, so goodness of fit of the initial wave is best gauged by the fit to the latter part of the data. this paper, we examine two questions: First, given its neglect of heterogeneities in population contacts, why does the well-mixed SIR differential equation fit real pandemic data so well? Second, when we fit the SIR model to real observational data, can we interpret its fitted parameters, when the interaction processes that are parameterized are themselves not realistic? We find that the fitted β parameter implicitly accounts for heterogeneous interaction dynamics, and the N tot parameter reveals an effective epidemic population size. We then apply our findings to the surprisingly well-fit omicron variant surge.

A. A Heterogeneous Extension
To examine why the simple SIR model fits real pandemic data at all, we consider a multi-population SIR model that allows for populations to be heterogeneous in their contacts and susceptibilities. We then interrogate simulations resulting from this model to examine how important such heterogeneity is in the dynamics, and what fitting such a model with a homogeneous model means mathematically.
The model assumes that the population is split into K subpopulations, with population sizes N i , i = 1, 2, . . . , K, that interact with each other with contact rate β ij , i.e. β ij is the number of infections a single infected individual in subpopulation i could cause in subpopulation j, if j contains only susceptible individuals. Our subpopula-tions represent geographic or demographic partitions of a population of size N tot . In this first treatment we neglect differences in recovery rate based on individual characteristics. The susceptible, infectious, and removed quantities of subpopulation i is given by S i , I i , and R i respectively and also write i S i = S, i I i = I, i R i = R, and i N i = N tot . Consider the model: To emphasize the similarity between the multipopulation and the single-population SIR model we write the equation in a form with a simple SIR part, with arbitrary coefficientβ plus a residual. We obtain a type of best SIR fit by minimizing the L 2 norm of the ratio of the heterogeneous residual to the homogeneous term on an arbitrary interval (t 1 , t 2 ), i.e. we minimize the size of the function This yields the unique minimum We explore this quantity for a certain class of subpopulation interaction structure in Section (II B) and (II C).

B. Intermediate Mixed/Unmixed Models
We first consider the case where subpopulations either do not interact or else interact at identical rates: Here δ ij is the Kronecker delta, b a mixing parameter ranging from 0 to 1, and A ij is an adjacency matrix describing which subpopulations interact. Assume initially A ij ≡ 1. When b is 0 the populations do not mix, and infections spread within but not between subpopulations, and when b = 1 the subpopulations mix completely, effectively merging into a single homogeneous population by construction of Eq. (5). Scaling β ij by N j ensures that interactions between subpopulations are proportionate to their sizes. A homogeneous SIR model fits this model for b 0.2 (Fig (4)). Moreover, theβ gathered by Eq. (4) is very Shown are the computed clusters grouped together by the methods described in Part I. The least-fitting cluster (in the authors' opinion) lining up with Arkansas still fits reasonably well when the fit at later stages is considered, especially when one keeps in mind that case data is highly unreliable and sensitive to stochastic effects at the beginning of the pandemic. close to the β achieved by least squares model fitting to the full simulation, although unsurprisingly least squares fitting performs better at fitting the actual simulated curve for lower values of b. The agreement betweenβ and the fitted β is encouraging: for the consideration of real data, one can only deduce model parameters via some sort of fitting algorithm without knowledge of subcompartmental dynamics, but the result suggests that least squares fitting optimally estimates the susceptibility parameter to maximize the simple SIR part of the model relative to the interpopulation dynamics.

C. A Model with Network Structure
To investigate possible effects of subcommunity interaction structure on the homogeneous-like dynamics, we now model subpopulations whose interactions are prescribed by a network with adjacency matrix A ij . A ij = 1 indicates two subpopulations that maintain frequent contact with one another, such as a pair of communities which go to the same grocery store or school. For the purpose of analysis, we modeled random connections between subpopulations as random Erdős-Réyni networks parameterized by mean degree. The SIR model approximates the graph dynamics model above b = 0.2, with the fit improving as the mean degree of each node is increased (Fig (5)). Just as for a complete graph,β agrees with β from least squares fits. Narrowing of error bars as mean degree of the random network increases express the decreasing importance of network structure as the network becomes more densely connected.
However,β does not asymptote to β = 0.2 whether we use Eq. (4) or least squares fitting. Instead, increas-ing mixing allows the model to gain awareness of how sparse the connections between subpopulations are, ultimately causingβ to decrease. Fig (5) suggests that as the mean degree increases, we expect that rate of decrease to slow, and thatβ will eventually asymptote to 0.2 as the network becomes complete (Fig (4)). Indeed, when all subpopulations have the same size,β may be written as Both numerator and denominator include a term which relies on graph structure and gains weight with the mixing parameter b and the covariance between the fraction of susceptibles and infecteds across subpopulations. The Si S ( Ii I − 1 K ) term represents a covari-ance measurement of the susceptible and infected percentage among the subpopulations, and is expected to be negative since, heuristically speaking, an increase in the number of infecteds corresponds to a decrease in the number of susceptibles, at least during a surge in cases.

The
term encodesβ's dependence on subpopulation interaction structure. For a complete network, this term is 1, and it decreases to 0 as the network becomes more sparse.
In the complete graph case and the case with network structure, we observe thatβ is dampened by heterogeneous community contacts. Givenβ's correspondence with the β gleaned from data fitting, we now assume that the SIR parameter fit for β actually underestimate the true person-to-person contact rate.

D. Comparison with other estimators of disease spread
We have shown that one can generate a single population-level transmission rate, β, for a heterogeneous population. From β we may derive the basic reproduction number R 0 = β/γ, the expected number of secondary cases produced in a completely susceptible population, by a typical infectious individual [23]. Our estimation method relies on minimizing the error between heterogeneous and an SIR model over the entire time course of a surge through Eq. (4). Prior estimation methods rely on the fact that in an SIR model, the parameter β describes the linearized, or initial, exponential rate of growth of the number of infectious individuals. For our model, we estimate this linearized rate of growth by two methods: 1. fitting an exponential on the first 20 days of the simulation, and 2. the next-generation matrix method from [24] (Fig (6)).
In [24], R 0 is computed for a general compartment model from the Jacobian matrix of the system. This matrix is evaluated at a disease-free equilibrium to determine the average number of individuals that a typical R0 estimates from fitting entire case data curve optimally to well-mixed model (blue) agree with linearized analysis by next-generation matrix method (yellow) [24], but not to empirical fits to the data assuming exponential growth (red).
infectious person infects when the population is asymptotically disease free, i.e. S i = N i for each i. Applying the method in [24] to Eq. (5) and making the assumption that N i = N j for all i and j as in the simulations, we obtain: where ρ denotes the spectral radius of the matrix A ij . In the complete network case (as in Fig (4)), ρ(A ij ) = K − 1, which yields R 0 = β/γ, matching the asymptote in Fig (6).
As the mean degree of our random network increases, R 0 values computed from the next generation matrix method in [24] and the R 0 computed fromβ converge (Fig 6). On all levels of network connectedness, and mixing parameter, exponential fitting consistently underestimates the contact rate. Under-estimation results from initial slowing in early growth, due to transmission being slower between subpopulations than within them. Conversely, for small values of the mixing parameter b, the next generation matrix method R 0 exceeds the optimal estimate fromβ, but the two estimators converge consistently at b values between 0.1 and 0.45, depending on the mean degree of connectedness between subcompartments.
The value of R 0 computed in [24] comes from the linearized dynamics; for example it provides a threshold for the stability of disease-free equilibria (see [25]). By contrast the estimate for R 0 computed in this paper is computed by approximating the spread of the disease by a homogeneous model. Surprisingly, the two methods produce confluent results even under modest levels of mixing between subpopulations.
In contrast to constant-parameter growth fitting, time series fitting calculates the time-varying rate of exponential growth of the number of COVID cases, and thence infers the number of new infections caused by each COVID case. Unlike the SIR model, R 0 , now called simply R, is not a constant, but typically varies over the course of a surge, and reflects not just the linearized dynamics of disease spread if an infected individual were transplanted to a population containing only susceptibles, but an actual estimate of new infections. Both SIR model and data fitting produce case curves that agree well with the surge of Omicron-variant cases, for which we use California as a representative example (section II F). However, the SIR model achieves this fit by assuming piecewise constant β, with a single (β, N tot ) pair covering most of the surge. We show in the next section that it is possible to forecast the end of the surge after an inflection point in a surge, and the total number of cases it will cause. Under the model, the surge ends only when the disease has been transmitted through an entire well-mixed subpopulation of size N tot . By contrast, the time-varying R value inferred by Epiforecasts [26] gives a compelling visualization of how transmission rates decrease during the surge. However, data fitting can not distinguish between decreases in transmission rate due to inevitable decrease of number of susceptibles around each infectious individual, or due to public health orders changing the course of the epidemic. The success of the SIR model in fitting the data weights the first factor over the second, though does not discount the effectiveness of public health measures implemented at or before the start of the surge in controlling its trajectory.

E. The Ntot Parameter as Effective Case Surge Size
In the derivation of the SIR model, N tot is the size of the population through which the disease is being transmitted. When using the model for data fitting, N tot is often treated as the size of the population; e.g. state or country from which the data was sourced [6]. However, real COVID cases occur in hot spots, and may not involve every individual in the studied population. In terms of our representation of this studied population by linked compartments, linkages between some compartments may be so weak that cases in one do not lead to a number of cases in the second that does not scale with the second compartment size. For this reason, we take advantage of the flexibility within our model of allowing N tot be fit alongside β.
Much like how we must reconsider the fitted β in the context of heterogeneities, we must now reconsider what the fitted N tot represents. The SIR model believes that N tot gives the scale of the population modeled. In turn, our fitted N tot parameter tells us something new and useful: the epidemic population scale that the data is conveying. Such a fitted N tot can then be a new, data-driven scale by which we can estimate per-capita case numbers in the context of a given surge.
Fitting N tot , however, is surprisingly tricky. In the SIR model, N tot always appears coupled in a fraction with β, and as such, when the scale is not necessarily apparent from the data, only the ratio β Ntot can be tuned to fit. The scale of the initial conditions used to fit the data are also unreliable, in part due to a large degree of uncertainty of what the initial conditions are given cumulative case data, i.e. we have little information about how many people can be considered susceptible or removed at any given time in order to seed the fit with an accurate initial condition. Indeed, when not enough data to determine scale is present, the fits for β and N tot are highly sensitive to the initial parameter estimates in our fitting algorithm (Fig 7). 7. Fitting analysis of the Ntot parameter for California's cumulative COVID case data during the Omicron variant surge. Left: the SIR model's best prediction (orange part of curve) when given a four different portions of the entire surge data (blue part of curve). Vertical green lines indicate the point where the SIR model starts predicting. As more data is given, the model better matches the actual cumulative case data (translucent black line). Right: the distribution of fitted Ntot parameter values as a function of number of days of data given (horizontal axis in plot). Fit distributions are generated by taking the negative exponential of mean squared residual of the fit after randomizing the initial parameter guesses (see section SID). The distribution narrows around a more unique guess as more data is given.
We analyzed omicron surges (which are fit quite well by the SIR model, as explored in the next section) to determine exactly how much data is needed for scale to be fit robustly (Fig 7). We observe that predictive power and parameter fit robustness rapidly increase immediately after the inflection point in surge data, suggesting that the model needs to be sure of a future plateau of cumulative cases to truly determine scale independently of the ratio β Ntot . This observation is supported mathematically by an analytical result on the eventual number of susceptibles as a function of model parameters: We see here a second indicator of scale for the SIR model present in the terminal dynamics of the model's trajectories. We reason that the inflection point in a case surge is the point after which these terminal dynamics become determined. In turn, the N tot parameter reveals itself uniquely.

F. SIR Curve Fits and Parameter Conclusions from U.S. State Omicron Case Surge Data
Armed with the more nuanced interpretations of what fitted β a N tot mean, we sought to draw conclusions from U.S. state data from the particularly infectious Omicron variant as it first rose to dominance among infectious individuals. Despite the relatively wide coverage of vaccines, the Omicron variant infected individuals as if they were purely susceptible, although vaccinated individuals were far less likely to end up hospitalized [5]. This high infection rate paired with the timing of two major U.S. holidays (Thanksgiving and Christmas) almost invariably led to homogeneous SIR-like dynamics for entire states (Fig 3), despite vast differences in population density and demographic structure. In fact, the only states not to exhibit clear homogeneous SIR dynamics appear to have extremely low-quality data. Given our early interpretations utilizingβ (Eq. 4), we conclude that some level of subpopulation interaction is happening state by state, with well-mixed dynamics occurring in the different subpopulations. Since the data considered is a surge of cases, we also expect the fitted β to represent the bulk susceptibility parameter dampened by subpopulation interaction effects due to negative covariance between susceptibility numbers and infectious numbers (Eq. 6). Our parameter fits for β in the following, therefore, are lower bounds for the true person-to-person susceptibility rate rather than actual estimates.
We chose to use California as a case study for our parameter fits due to the absolute certainty of there being at least two loosely connected population centers for COVID spread (the San Francisco area and the Los Angeles area, separated by approximately 400 miles). Much like the vast majority of states, California's Omicron surge is extremely well approximated by a simple SIR curve. The most frequently-occurring set of parameters that fit the data put the basic reproduction number at approximately R 0 ≈ 3 (a lower bound just as β is), and N tot ≈ 12 million. The SIR model therefore thinks that California's data is actually represented by homogeneous spread among 12 million individuals. Our earlier interpretations temper this conclusion with the reality that the spread cannot be homogeneous, and leads us to the more nuanced conclusion: the spread of Omicron in California is dominated by the well-mixed disease spread amongst 1 or more distinct subpopulations totaling to approximately 12 million individuals, around 30% of California's population. Among these communities, spread has a very rapid rate of at least R 0 = 3. This highlights the dominant role that a closely interacting set of individuals can play in driving pandemic case trajectories.

III. DISCUSSION
In the Part I of this paper we showed that COVID case growth curves from different US states can be collapsed onto a small number of universal curves, after translation in time and rescaling of the total number of cases. These symmetries are present in the SIR model, and indeed we find, alongside previous works [6,27] that SIR models can well fit the initial COVID surges in 26 of 52 states and territories in our data set, although our SIR model fits use the simplest, constant-parameter ODE version.
Our analysis of data from states and territories provide empirical evidence that there is enough mixing between the disease hotspots within each state or territory that approximation by a well-mixed model is appropriate. The concordance between β obtained by curve fitting and Eq. (4) affirms that least squares fitting extracts information on the subpopulation interaction structure and dynamical asymmetry between subpopulations, and highlights the SIR model as a coarse-grained model for disease transmission in heterogeneous populations.
Remarkably, under moderate levels of transmission between subpopulations, our model-fitting approach, which is based on fitting the entire disease case curve, produces parameter estimates that agree well with the next generation matrix model based only on the linearized dynamics that include information only on initial epidemic features. We additionally gain the ability to reinterpret the total epidemic size N tot as a new, data-driven normalization factor.
We may therefore construct case rates with respect to surge size, effectively considering normalized cases only among the population actively transmitting.
In general, simplicity of fitting to real data makes SIR (and similar well-mixed models, such as SEIR) powerful tools for predicting the ongoing course of an epidemic. However, the assumptions behind the model drastically simplify real patterns of human interaction and disease transmission, and the question of estimating community interaction structure remains. Indeed, an important corollary of our analysis is that, with even modest levels of connectedness between the subpopulations, they can function as well-mixed, making the details of the substructure undetectable by a well-mixed model. Although heterogeneity-capturing models aim to render these interactions, they introduce additional parameters that often must also be fit to the data. It can be hard to distinguish improvements in fitting due to greater model realism from improvements due to increased parametric flexibility.
Data fitted analyses make no assumptions about the underlying patterns of interaction, but create short term predictions based on fitting evolving exponential growth curves to the data. SIR models do not fit every phase of the case data; for example we could not fit the second surge identified in Part I, likely because of the complicating effects of time evolving contact rates caused by chang-ing work patterns and public health orders, and the presence of previously recovered individuals within the population. By contrast, although models of the the Omicron surge encounter similar complications, reluctance to reimpose social distancing measures, and the ability of Omicron to readily infect vaccinated and previouslyinfected individuals lead to a time course that is closer to the first surge, and that can be well-fit by an SIR model.