We developed and validated single and multiple population ordinary differential equations (ODE) and agent-based models of COVID-19 spread within and between defined groups of individuals. The general model logic (Multimedia Appendix 1) was adapted from classic susceptible-exposed-infected-recovered (SEIR) frameworks [67,68] similar to those described elsewhere for many other infectious diseases [63-76]. The single population model describes spread dynamics within one defined population (eg, on-campus students or local community residents), whereas the multipopulation model additionally includes cross-exposure between two or more groups. Multiple change points were included for all parameters to allow for policy or behavior changes when fitting models to historical data.

State variables at time t include the numbers of individuals in population j that are free of and susceptible to COVID-19 (S_j(t)), exposed to COVID-19 but not yet infectious themselves (E_j(t)), COVID-19–positive and infectious to others (I_j(t)), recovered and not susceptible to reinfection (R_j(t)), and COVID-19–associated deceased (D_j(t)). Exposed individuals are assumed to have a small chance (p) of warding off an infection before becoming infectious. Recovered individuals are assumed not able to be reinfected within at least a 16-week (one semester) time frame [37,77]. Each state variable is updated numerically at each time increment (set here to 0.01 days) based on its previous value, values of other state variables at the previous time step, and the equations governing their interdependent relationships, with this process continuing iteratively for 16 weeks.

For example, the number of individuals in the susceptible population (S_j(t)) is decremented by the number of newly exposed individuals (S_j(t) · expo_j) and increased by the number who previously were exposed but did not develop infections, (E_j(t) · expr_j), where the daily exposure rate expo_j, the average risk of transmission multiplied by the average number of contacts per day, is back-computed from the basic reproduction number R₀ (average number of new infections per infected individual) and recovery and mortality rates, and the recovery rate of noninfectious exposed individuals expr_j is the inverse of the corresponding recovery time t_{rec, nonj}.

In turn, the number of exposed individuals is increased by S_j(t) · expo_j and decremented by the number who develop infections (E_j(t) · infe_j), where the daily infection rate infe_j is the ratio of the probability of becoming infected upon exposure p_jover the average incubation time t_incj. Infected individuals either recover or die at rates of infr_j · I_j(t) and mort_j · I_j(t), respectively, where the daily recovery and mortality rates are the inverse of the average recovery time t_{rec, incj} and the ratio of the overall COVID-19 case fatality rate for that population (CFR_j) over the average time from infection until death t_i2dj, respectively.

The governing rate change dynamics for each state variable at each time step during numeric evaluation thus are the following:

S_j(t) (susceptible; + not-infected/infectious (nor immune) after exposure – new exposures due to within-population and between-population contact with infectious individuals):

E_j(t) (exposed; + new exposures – past exposures now infected/infectious – past exposures now not infected/infectious [now susceptible]):

I_j(t) (infectious; + past exposures now infected/infectious – past exposures now not infected – deaths):

R_j(t) (recovered; + infected individuals who recover [with immunity]):

D_j(t) (deceased; + COVID-19–related deaths):

where

(rate at which people transition from susceptible to exposed)

(rate at which people transition from exposed to infected)

(rate at which people transition from infected to recovered)

(rate at which people transition from exposed to recovered)

(rate at which people transition from infected to deceased) and

(sum of all subpopulations in region i at time j)

where χi,j = 1 if populations i and j interact and 0 otherwise and p_j = the proportion of exposed individuals that transition to infected (versus recovering to susceptible). The multipopulation models allow for separate parameter values for each population, such as based on their demographics, with a cross-exposure parameter (ri_j) defining the relative rate at which infectious individuals in one population expose susceptible individuals in the other (typically lower than within-population, assuming less interaction).

Model accuracy was validated using standard methods [65,78-82], cross-validation, and varied state and county empirical data (January to July 2020) exhibiting different epidemic patterns, magnitudes, and timings (Multimedia Appendix 2). Model results closely emulate historical data across multiple settings, with accuracy on par with or exceeding norms reported elsewhere [83-89] and with ≤1 change points generally providing good fits, suggesting good prospective short-term prediction capability.

Model inputs (Table 1) used in the community and campus models were estimated using a combination of published and grey literature, expert opinion, and search-based optimization. For campus inputs with uncertainty, we used Monte Carlo simulations to create 1000 synthetic results across plausible ranges, using the shown most likely, maximum, and minimum values to generate asymmetric triangular distribution random variates. Since little data exist about on-campus spread [21], for exposure rates we used the shown ranges for the average number of infected students divided by the exposed-to-infectious percentage.

For community populations, we further calibrated inputs via a particle swarm search algorithm to minimize root mean square error differences between historical and model-predicted infections and mortality, running each parameter search 1000 times. For model fits with change points, separate values for all inputs were optimized for each time segment, with state variables at the start of each new time segment set to their values at the end of the prior segment.

For initial disease prevalence in the local community and among arriving students, we also used expected values and probability intervals from a logistic growth curve fit to historical COVID-19 infection counts, estimating prevalence among arriving students at the start of the fall 2020 semester using a weighted average of prevalence predictions based on home locations. Resulting community and student prevalence ranges were validated against data reported in the media. Positive individuals at the start of the semester were assumed distributed between exposed but not yet infectious (24.3%) and infectious (75.7%) groups based on approximate relative durations that an average infected individual might spend in each state. All model inputs were based on published literature listed in Table 1 or aggregate state and county infection and mortality online data [103] and thus not subject to human subjects internal review board approval.

Common university reopening scenarios were identified from literature and published surveys [7,19], generally belonging to one of several categories (see Table 2 for examples), which then were used to estimate plausible ranges for R₀ reductions. The most common approaches included primarily or fully in-person (35%), primarily or fully online (32%), and hybrid (19%) [7]. As examples, the University of Washington reopening plan [104] exemplifies a conservative approach, with more than 90% of courses taught online, courses relying on direct interactions (eg, medical and health sciences) taught in person with safety precautions, the majority of student services and advising taking place remotely, and any staff who can work remotely doing so. In contrast, Purdue University illustrates an opposite approach [105], with classes mainly taught on campus with contact precautions until the Thanksgiving break, relying on students to manage their personal safety.

As most reopening plans involve reducing either interpersonal contact or infection spread, we implemented these as multiplicative reductions in the reproduction number R₀, with effect sizes of individual actions estimated via literature estimates and expert opinion. Overall scenarios then were defined with estimated effects on R₀ spanning the base case of no change, small 25% cumulative reductions (eg, Rice University, Ohio State), moderate 50% reductions (eg, Northeastern University, Connecticut State), and 75% reductions as a best-case scenario for comparison, with each scenario coupled with initial student infection rates of 0.1%, 1%, 2%, and 5% based on university reporting.

For each scenario, 1000 model replications were run for campus alone, community alone, and campus × community combined to estimate additional cross-exposure impacts of each population on the other. For the campus × community cases, each of the 1000 community parameterizations were randomly coupled with the 1000 random sets of campus inputs, with the two populations interacting via 1000 random values of the cross-exposure parameter, ri, sampled from the range shown in Table 1.

Given that the ratio of campus-to-community population sizes may affect cross-infection results and prevention policies, we assumed three general settings: (1) an urban campus of 10,000 students with 100,000 residents living in the immediately surrounding residential areas or neighborhoods (in which off-campus students tend to reside); (2) a student body of the same size (10,000), but with fewer (40,000) residents living close to the campus; and (3) a smaller number of 2000 students with 40,000 residents living near campus. The first scenario might represent a large university in a major city, whereas the second might represent a large rural university, and the third a smaller undergraduate college in a nonurban setting, although these student-to-community populations (1:10, 1:4, and 1:20) can be extrapolated to other settings with similar ratios.

Monthly and total counts of COVID-19 exposures, infections, and deaths for each population were tabulated and plotted longitudinally. First inflection points (dates of steepest increases) for each outcome, scenario, and population were identified numerically, since in diffusion theory interventions after these points tend to be less effective. To estimate cross-exposure effects, pairwise differences were computed between each of the 1000 campus results and their campus × community counterparts, and similarly between each of the 1000 community results and their campus × community counterparts. For all model results, medians, standard deviations, and 95% probability intervals were computed, with 10% trimming to reduce any extreme outlier replicate effects.

Sensitivity analyses were conducted to identify model inputs to which the mean and variance of results are most sensitive via central composite factorial experimental designs [116] as this could inform policy-making, interventions, and target setting. Results were analyzed using general linear models including linear and pairwise interaction terms for each outcome (replication means and variances of total and additional campus and community infections and deaths), with resulting effect coefficients normalized to their corresponding ranges and ranked according to statistical significance.

Within any given assumptions for COVID-19 prevalence among arriving students and semester initialization precautions, the predicted number of students per 10,000 who might be exposed, be infectious, and die over a 16-week semester could vary by up to 10-fold (Figure 1). By semester end, under the base case (2% arrival prevalence, little returning precautions and/or effectiveness) predicted student outcomes range from 471-9458 infections (median 2286, SD 2627) and 0-123 deaths (median 9, SD 14).

The more realistic case (1% arrival prevalence) reduces these consequences to a median of 1332 (SD 2552) infections and 5 (SD 12) deaths, with the steepest increases in exposures and infections typically occurring at midsemester onward (with important implications on geographic spread as students return to their home communities). Although less likely, idealistic (0.1%) and worst-case (5%) initial prevalence scenarios were also considered for comparison, given that epidemic prevalence might change in future semesters. The first would result in a median of 158 (SD 1760) infections and 1 death (SD 6 deaths) by semester end and the latter in 3996 (SD 2485) infections and 16 (SD 17) deaths.

Note that the left-hand plots in Figure 1 depict the current number of exposed people on any given day (with each exposure spanning several days) in order to give an indication of the changing amount of contract tracing and isolation required, as well as the changing population risk, whereas the other plots depict the cumulative number of infections and deaths to date in order to summarize the total public health impact. Under the base case scenario, the number of active student exposures at any given time (eg, for contact tracing and isolation) ranged from 3-1576 per 10,000 individuals (eg, as high as 15% of a student population), with significant implications on resource planning and viability (Table 3). Under the two most likely scenarios, by midsemester the total number of infections might be as high (mean plus one standard deviation) as 810 per 10,000 students (with 4 deaths) or as low (mean minus one standard deviation) as 782 per 10,000 students (with 0 deaths).

Similar results occur assuming the various semester precautions summarized in Table 2, corresponding to plausible reductions in R₀ of 25%-50% (Figure 2). In general, current strategies to reduce exposure during a semester appear effective, although under most scenarios a concerning number of students still can become infected or die. Even in the very optimistic case of a 75% R₀ reduction, included for comparison as a hypothetical “best case” scenario, 95-132 infections (median 107) and 1 death per 10,000 students may occur by midsemester, increasing to 97-139 infections and 0-3 deaths by semester end.

Figure 3 summarizes additional community (blue lines) and campus (red lines) impacts of reopening due to campus × community cross-exposure, assuming the same scenarios described above; for comparison, the top row shows the baseline number of community exposures, infections, and mortality without reopening. Local community impacts (Table 4) of opening with little-to-no semester operation precautions and/or adherence might range from 1-9768 additional community infections (median 158, SD 1131) and 0-491 additional community deaths (median 6, SD 53).

The two more realistic scenarios result in a total of (for 25% exposure reduction) 1-5577 additional community infections (median 56, SD 516) and 0-272 additional community deaths (median 3, SD 24), and (for 50% exposure reduction) 0-464 additional community infections (median 14, SD 45) and 0-23 additional community deaths (median 1, SD 2). For comparison, the hypothetical best-case scenario with 75% exposure reduction results in 0-33 additional community infections (median 2, SD 4) and 0-2 additional community deaths (median 0, SD 0.2).

The corresponding impact of the community on student outcomes is smaller, ranging from 0-390 additional student infections (median 21, SD 38) and 0-2 additional student deaths (median 0, SD 0.2) for the base case to 0-42 additional student infections (median 5, SD 6) and 0 additional student deaths (median 0, SD 0.04) for the idealistic case (75% exposure reduction). The two more likely cases result in (for 25% exposure reduction) 0-279 additional student infections (median 17, SD 30) and 0-1 additional student deaths (median 0, SD 0.1), and (for 50% exposure reduction) 0-115 additional student infections (median 8, SD 12) and 0-1 additional student deaths (median 0, SD 0.06).

To estimate the impact of school size and location, Figure 4 compares results under other student-to-community population sizes, assuming the same arrival prevalence and campus operation precautions, compliance, and effectiveness scenarios. While intuitive differences exist in raw totals, results are similar and scale-invariant after adjusting for population size. For example, multiplying results for the second case of 40,000 residents by 2.5 yields similar curves to those for the first case of 100,000 residents. This suggests that the above results may generalize to other settings and that between-location differences in epidemic patterns (and therefore in public policies to limit spread) likely arise from variations in campus × community interaction, rather than in population ratios.

Finally, factorial sensitivity analysis produced the relative parameter rankings shown in the two rightmost columns of Table 2, which follow intuition and serve as further model validation. The most statistically significant factors (main effects) affecting the expected (Multimedia Appendix 3) and variation (Multimedia Appendix 4) in total campus infections were recovery time, R_0,campus, incubation time, initial prevalence among arriving students, and R_0,community. Similarly, for total community infections, the most important factors were recovery time, R_0,community, incubation time, initial prevalence among community residents, and community case fatality rate. Numerous interaction terms were also significant in both cases, as would be expected in such a model.

For additional community infections, no dominant factors that affect either magnitude or variability of outcomes were evident beyond the basic reproduction numbers of each population (R_0,campus, R_0,community) and the campus × community cross-exposure rate ri. The large number of other statistically significant main effects and interaction terms also underscore the multidimensional challenge of predictably limiting community impact and by extension the importance of effective surveillance and mitigation.

The COVID-19 pandemic continues to be a significant public health crisis, with infections and mortality in many regions meeting or exceeding those in early 2020 before physical distancing and closures were implemented. With many colleges and universities reopening, model-based analyses can help inform these important decisions as well as the degree of uncertainty in the resulting outcomes. Three important results of the present analysis are the following: (1) infections and mortality from campus reopening are highly variable and nearly impossible to predict with any certainty, (2) reopening campuses can significantly impact local communities even under best-case scenarios, and (3), while few exist, prevention and public health measures that target campus × community exposure could be effective.

While conditions may exist under which reopening is relatively safe to the local community, at present these appear in the significant minority. Our results also agree in general with emerging empirical data from the fall 2020 semester, including reports that COVID-19 deaths in US communities with open colleges roughly doubled from August to December 2020, compared with a smaller 58% average increase in communities without colleges. Genetic sequencing results further suggest that many deaths in college towns were of older people who had contact with infected students [117].

Several important public health implications of our results exist. First, decisions about whether to open in future academic terms or epidemics should be informed by updated model inputs, projections of local conditions, and campus × community public health measures. Second, since any trajectory within the produced intervals could occur, reopening decisions should consider these ranges rather than averages alone. Third, given the wide uncertainties in results from reopening, criteria should be established for rapidly detecting when to tighten precautions. Fourth, contact tracing and isolation capabilities should be ensured to be sufficient to respond to the range of model results.

Like any model-based analysis, results herein have some limitations and simplifications. A common barrier in such models is data availability for input estimation and results validation (hence our search-based approach). The deterministic ODE modeling framework ignores inherent variability and population heterogeneity [118], motivating our use of Monte Carlo analysis, parameter search replicates, and randomly sampled scenarios. Standard model simplifications include limiting the number of populations (eg, one overall homogenous community population), limiting spread to just SARS-CoV-2 (eg, ignoring seasonal influenza, substance abuse [75,76], and co-epidemic impacts), and not time-varying precaution compliance as concerns and vigilance relax or heighten over time. Some scenarios were also included for potential insights rather than being feasible in practice (eg, 75% reduction in R₀, near 100% precaution compliance).

Further work could expand on these results, including addressing some of the above simplifications, rerunning analyses for future semesters using more recent data for model calibration, and considering more heterogeneity in community and student populations. Future work could also seek to determine combined conditions (reduced prevalence, vaccine effectiveness, improved precaution methods, etc) under which outcomes are both safer and more certain. Public health reopening and precaution decisions at citywide or statewide levels also might be examined, such as alternating on-campus semesters or limiting combined student densities, to manage net community risks.

Controlling the COVID-19 pandemic is extremely critical. Mathematical models can offer valuable insights to inform important public health and policy decisions, including potential community and campus impacts from university reopening. The analysis summarized herein suggests that outcomes over a 16-week semester can be highly unpredictable under any set of assumptions or precautions, with three important implications: (1) community impacts from campus reopening are highly difficult to predict in advance, (2) on- and off-campus surveillance and response methods therefore are critical, and (3) additional precautions to reduce impacts of open campuses on local communities appear warranted.