1 Introduction

In late 2019, a novel coronavirus, SARS-CoV-2, was identified as the cause of a cluster of pneumonia cases (Zheng-Li 2021). On March 11, 2020, the World Health Organization declared the 2019 novel coronavirus outbreak (COVID-19) a pandemic (World Health Organization 2020). Shortly after, nearly every higher-education institute rapidly transitioned all classes to on-line instruction to “flatten the epidemic curve.” As of October 13, 2022, the cumulative number of confirmed COVID-19 cases exceeds 620 million (Johns Hopkins Corona Virus Resource Center 2022). Although the availability of multiple effective vaccines offers the likelihood of a return to normal life, the emergence of highly infectious variants and the advent of booster shots means that the return to our pre-COVID-19 existence is not in our immediate future (Mathieu et al. 2021).

Starting in Fall 2020 in the USA, colleges and universities have attempted to employ strategies to manage COVID-19 with mixed results. Overall, there were substantial increases in the number of new COVID-19 cases after school re-opening (Nierenberg and Pasick 2020). Moreover, even though, by age, college students are less likely to experience severe complications from COVID-19, the same is not true for their surrounding communities. During the winter of 2020, large surges in COVID-19 cases from college students were followed by subsequent infections and deaths in the wider community (Ivory et al. 2021). In addition, many campuses delayed their in-person instruction in early 2022 due to emergence of the omicron variant (Moody 2022). While there is a strong desire for higher-educational institutions to maintain in-person instruction, it is clear that for the foreseeable future this will require an effective COVID-19 management policy.

Nationally, educational institutions need to evaluate how to most effectively plan activities during the academic year while ensuring they do not contribute to local outbreaks (Ellis 2020; Paltiel et al. 2020; Wrighton and Lawrence 2020). Although some campuses, such as the University of California (UC) and California State University systems, are mandating the COVID-19 vaccine for all students and employees, these mandates will not be required by all campuses or campus populations (Callimachi 2021). In places where the COVID-19 vaccination is not mandated, the population vaccination levels are likely to vary with local COVID-19 vaccine acceptance patterns (Lazarus et al. 2021).

Mathematical models have a proven track record of providing novel insights into the spread and control of epidemics. Dynamic epidemic models have been used to study COVID-19 at many scales (Chinazzi et al. 2020; Gatto et al. 2020; Gilbert et al. 2020; Hethcote 2000; Kucharski et al. 2020; Lofgren et al. 2021; Mizumoto and Chowell 2020). Given the wide-spread campus closures due to COVID-19, models have been developed to study the spread of COVID-19 on college campuses to evaluate reopening strategies (Gressman and Peck 2020; Lopman et al. 2020; Weeden and Cornwell 2020). Many such models have considered only student populations, thereby reducing the population to an isolated “bubble” in which students do not interact with individuals outside the campus and faculty are not accounted for. In this study, we develop a model for COVID-19 dynamics in “bubble-like” institutions—such as universities, nursing homes and prisons—where individuals within those communities have complex structured interactions defined by their roles, but, rather than a bubble, the boundaries between these environments are porous and certain types of individuals—professors, staff, guards—intermix freely within a larger surrounding community where COVID-19 is also spreading. We began the modeling work we present during the summer of 2020, as our campus and other colleges around the world worked to manage instruction. While our initial underlying framework was developed before the emergence of multiple variants, we have modified our model to include the possibility of reinfection for both vaccinated and unvaccinated subpopulations.

In this study, we develop a structured SEIR model of COVID-19 dynamics on a college campus and investigate the sensitivity of behavior to the vaccinated population on campus and other non-pharmaceutical interventions (NPIs) such as mask-use and social distancing. Our goal is to understand how vaccine hesitancy both within the campus population and the surrounding community will impact disease propagation and which interventions will be the most effective. More specifically, we individually model the various subpopulations at the university, including on-campus undergraduates, off-campus undergraduates, graduate students, and faculty/staff. We connect our campus to the surrounding community where behavior outside the university will impact COVID-19 dynamics within the university. We perform a global sensitivity analysis of model behavior—cumulative number of infections at the end of the semester and infection doubling time—and consider the first-and total-order effect of epidemic parameters and social contact behavior.

In Sect. 2, we first develop our structured epidemic model, then describe the model outputs we will study as well as the variance-based sensitivity analysis approach we employ. In Sect. 3, we discuss the campus data we use to parameterize our model. Although we use the campus network of UC Merced, we believe our results are representative of mid-size rural colleges. In Sect. 4, we present our results. We conclude in Sect. 5. We note that under all conditions, NPIs are still important to mitigating the spread of COVID-19 on campus and urge universities to continue to support their use.

Fig. 1
figure 1

(Color figure online) a Contacts between the 4 campus populations (on-campus undergraduates (d), off-campus undergraduates (u), graduate students (g), faculty/staff (f)) and the outside community. Contacts are separated into classroom (black), dormitory (black dotted), off-campus housing (green), and the outside community (orange). Thickness of arrows corresponds to the number of contact hours. b The stages of the COVID-19 infection, included in the ODE model, that an undergraduate living off-campus would progress through: susceptible (S), exposed (E), asymptomatically infectious (\(I^a\)) or symptomatically infectious (\(I^s\)), if symptomatically infectious individuals can choose to self-isolate (at home) (H) or not (N), and finally, both asymptomatically and symptomatically infectious individuals recover (R). Note that some percentage of the population is initially vaccinated (V) and, for simplicity, we assume no more individuals will get vaccinated once the semester begins. Both vaccinated individuals and recovered individuals can become susceptible after certain period of time

2 Methods

2.1 Model Description

In this section, we describe the ordinary differential equation (ODE) compartment model that we have used to make most of the analysis for the university as presented to the administration during the summer and fall of 2020, see Fig. 1a. This SEIR model is unique as we model the four subpopulations: (1) undergraduate students living off campus, u, (2) undergraduate students living on-campus and in dorms, d, (3) graduate students, g, and (4) faculty and staff, f. Individuals are designated by both population and COVID-19 status. Figure 1b presents the phases of the diseases for one of the subpopulations, the undergraduate students living off-campus (u).

Each subpopulation in our model is divided into eight compartments related to different stages of infection and immune status: susceptible, S; exposed, E; asymptomatically infectious, \(I^a\); symptomatically infectious, \(I^s\); symptomatically infectious and in self-isolation (at home), H; symptomatically infectious but not in self-isolation, N; recovered, R; and vaccinated, V. A susceptible (S) individual may become exposed (E) to SARS-CoV-2 (the virus that causes COVID-19) by coming into contact with any infectious individuals (I) in any of the 4 subpopulations or the outside community (off campus). All exposed individuals become infectious, either asymptomatically (\(I^a\)) with probability \(\phi \) or symptomatically (\(I^s\)) with probability \(1-\phi \), after spending an average of \(1/\sigma \) days in exposed state. After an average of \(1/\gamma ^s\) days, symptomatic individuals may choose to self-isolate and/or report to health services for testing (H) with probability \(\alpha \) or not (N) with probability \(1-\alpha \). Finally, as we do not model disease mortality, all symptomatically infectious individuals eventually recover (R)—assuming the recovery rates for H and N individuals are h and \(\delta \), respectively. For asymptomatic individuals, we assume they recover after \(1/\gamma ^a\) days. For simplicity, we assume that both recovered individuals (R) and vaccinated individuals (V) are unable to contract the virus nor pass it to susceptible individuals for a certain period of time due to naturally acquired immunity (\(1/\kappa ^R\)) and vaccine-mediated immunity (\(1/\kappa ^V\)), respectively. We assume that the asymptomatically infectious individuals have lower probability in transmitting the disease than symptomatically infectious individuals, thus we employ the fraction of \(\beta \) of asymptomatically infectious individuals \(\zeta <1\) where \(\beta \) denotes the transmission rate for symptomatically infectious individuals. Figure 11 in Appendix A shows all the subpopulations and all the phases of the disease for each subpopulation.

For a subpopulation \(i\in \{u,d,g,f\}\), the dynamics of each stage of the infection are governed by

$$\begin{aligned} \frac{\hbox {d}S_i}{\hbox {d}t}&= -m \beta S_i F_i + \kappa ^V V_i + \kappa ^R R_i, \end{aligned}$$
(1a)
$$\begin{aligned} \frac{\hbox {d}E_i}{\hbox {d}t}&= m \beta S_i F_i - \sigma E_i, \end{aligned}$$
(1b)
$$\begin{aligned} \frac{\hbox {d}I^a_i}{\hbox {d}t}&= \phi \sigma E_i - \gamma ^a I^a_i, \end{aligned}$$
(1c)
$$\begin{aligned} \frac{\hbox {d}I^s_i}{\hbox {d}t}&= (1-\phi ) \sigma E_i - \gamma ^s I^s_i, \end{aligned}$$
(1d)
$$\begin{aligned} \frac{\hbox {d}H_i}{\hbox {d}t}&= \alpha \gamma ^s I^s_i - h H_i, \end{aligned}$$
(1e)
$$\begin{aligned} \frac{\hbox {d}N_i}{\hbox {d}t}&= (1 - \alpha ) \gamma ^s I^s_i - \delta N_i, \end{aligned}$$
(1f)
$$\begin{aligned} \frac{\hbox {d}R_i}{\hbox {d}t}&= \gamma ^a I^a_i + h H_i + \delta N_i -\kappa ^R R_i, \end{aligned}$$
(1g)
$$\begin{aligned} \frac{\hbox {d}V_i}{\hbox {d}t}&= -\kappa ^V V_i, \end{aligned}$$
(1h)

with initial conditions

$$\begin{aligned} \begin{array}{l} I_i^a(0)= I_{i,0}^{a} , \quad I_i^s(0)= I_{i,0}^{s}, \quad V_i(0) = \left( n_i -I_{i,0}^{a}-I_{i,0}^{s} \right) v_{i,0}, \\ S_i(0)= n_i -I_{i,0}^{a}-I_{i,0}^{s} -V_i(0), \\ E_i(0)=H_i(0)= N_i(0)= R_i(0)=0. \end{array} \end{aligned}$$

Here \(n_i\) denotes the total number of individuals in subpopulation i, and \(v_{i,0}\) denotes the percentage of subpopulation i that has been vaccinated at time \(t=0\). Note that all the subpopulations are coupled through the force of infection, \(F_i\), where an individual in any subpopulation can be infected by any infectious individual in the campus population or through an infectious member of the community. The calculation involves the contact matrix \(\mathbb {C}\),

$$\begin{aligned} \begin{bmatrix} F_u\\ F_d\\ F_g\\ F_f \end{bmatrix} = \mathbb {C} \times \begin{bmatrix} (\zeta I^a_u + I^s_u + N_u)/n_u\\ (\zeta I^a_d + I^s_d + N_d)/n_d\\ (\zeta I^a_g + I^s_g + N_g)/n_g\\ (\zeta I^a_f + I^s_f + N_f)/n_f\\ \zeta \psi M + M(1 - \psi ) \end{bmatrix}. \end{aligned}$$
(2)

The contact matrix \(\mathbb {C}\) is a \(4\times 5\) matrix, with the indices 1, 2, 3, 4, and 5 correspond to off-campus undergraduate students, on-campus undergraduate students, graduate students, faculty/staff, and outside community, respectively. Taking off-campus undergraduate students as an example, we can obtain the explicit form of the force of infection \(F_u\) from Eq. (2):

$$\begin{aligned} F_u= & {} c_{1,1} \dfrac{\zeta I^a_u + I^s_u + N_u}{n_u} + c_{1,2} \dfrac{\zeta I^a_d + I^s_d + N_d}{n_d} + c_{1,3} \dfrac{\zeta I^a_g + I^s_g + N_g}{n_g} \\{} & {} + c_{1,4} \dfrac{\zeta I^a_f + I^s_f + N_f}{n_f} + c_{1,5} (\zeta \psi M + M(1 - \psi ) ), \end{aligned}$$

where \(c_{i,j}\) denotes the (ij)-th entry of the contact matrix \(\mathbb {C}\). The first four terms signify the force of infection from off-campus undergraduate students, on-campus undergraduate students, graduate students, faculty/staff, respectively, and are generated in the same fashion. For the first term, we can break it into the product between the number of contact-hours the off-campus undergraduate student has with other off-campus undergraduate students per day \((c_{1,1})\), the probability of disease transmission per contact-hour (e.g. \(\zeta \) if asymptomatic, 1 if symptomatic), and the proportion of contacts that are infectious (e.g., \(I^a_u/n_u\) if asymptomatic, \((I^s_u+N_u)/n_u\) if symptomatic). The last term is generated by contact-hours between the average undergraduate student and the outside community per day (\(c_{1,5}\)) multiplied by the static population of infectious individuals in the surrounding community (both asymptomatic and symptomatic), where M denotes the probability that an outside community individual is infectious and \(\psi \) denote the probability that an outside community individual with COVID-19 is asymptomatic.

Table 1 Description of parameters and their ranges used in the global sensitivity analysis. See Appendix A for details on parameterization

The model has two types of parameters: (1) parameters related to COVID-19 epidemiology and (2) parameters related to contact patterns between the 4 subpopulations and the outside community (contact matrix \(\mathbb {C}\) is described in Table 2 and explained in Sect. 3). Details of the first type of parameters are given in Table 1. Details of contacts are given below in Sect. 3. The critical term in our model is the force of infection, defined in Eq. (2), which governs the spread of COVID-19 and is impacted by interventions such as mask-use, quarantine of infectious individuals, and changes in social connectivity, including class-size and housing caps through both types of parameters. We account for the effect of masks with our mask efficiency parameter which, when used in classrooms, reduced the transmission rate \(\beta \) by a factor \((1-m)\). In other words, \(m=1\) or \(1-m=0\) in Eqs. (1a, 1b) corresponds to our baseline case when masks are not required or wearing a mask has no effect on the transmission rate \(\beta \). We account for a quarantine period (1/h) and the probability that a symptomatic individual chose to self-isolate (\(\alpha \)). Finally, changes in social connectivity, such as changes in class-sizes, are implemented by modifying the appropriate contact matrix.

Because our goal is to investigate the sensitivity of COVID-19 dynamics to NPIs in the context of a vaccinated campus population, we consider two simplifying assumptions allowing us to include vaccinated individuals in each of our campus subpopulations. First, our goal is to model COVID-19 in the context of a university environment. As such, we followed the requirements of the UC, which requires vaccination before the semester begins and did not consider an on-going vaccination program. Thus, in our model we include all vaccinated individuals in the vaccinated compartment at time \(t=0\), and there is no inflow to the vaccinated compartment throughout the simulation. Second, we assume individuals would gain temporary immunity from vaccination which wanes over time. For simplicity, in our model, vaccinated individuals move to the susceptible compartment at the rate \(\kappa ^V\). Similarly, we assume the immunity recovered individuals acquired from infection also wanes over time and thus, in our model recovered individuals transition to the susceptible compartment at the rate \(\kappa ^R\).

2.2 Infection Doubling Time

We are interested in analyzing the effect of intervention strategies and vaccination on the resulting dynamics of COVID-19 infections on our campus. In our work, the infection doubling time, \(\Delta T\), is the characteristic number of days for the cumulative number of COVID-19 infections to double, \(C(t_{i+1})=2C(t_i)\), where \(t_{i+1}=t_i+\Delta T\) for \(i\ge 0\) and \(t_0\) is a point in time at the beginning of the semester, see Fig. 2. This quantity is a characteristic of the disease dynamics during the early stages of disease spread (beginning of the semester), when the number of on-campus infections remains low, and the susceptible population remains large (Patel and Patel 2020). This characteristic captures the start of a semester when students that are allowed back to campus have been tested for COVID-19 (large susceptible population), and the probability that an infectious student returns to campus is low. Infection control measures aimed at “flattening the curve” means increasing the infection doubling time (Nunes-Vaz 2020). In this work, we use the infection doubling time as a measure of epidemic dynamics to asses which intervention strategies are associated with increased variance in the infection doubling time.

Fig. 2
figure 2

(Color figure online) Cumulative Infections and Infection Doubling Times. The figure illustrates the doubling of the cumulative number of infections with respect to the infection doubling time (\(\Delta T\)). Here, \(C(t_i)\) is the cumulative number of infections at a time \(t_i\) since the beginning of the semester

During this initial period of disease spread, at the beginning of a semester, we observe an exponential growth phase in the number of infections (Fig. 2). If the number of infections is growing at a rate r, the infection doubling time is given by \(\Delta T=\log (2)/r\) (Muniz-Rodriguez et al. 2020). In this early stage of approximately exponential growth in the number of cumulative infections, C(t), where \(C'(t)\approx r C(t)\), we can estimate r by using cumulative infections at two points in time \(C(t_1)\) and \(C(t_2)\), where \(r=\log \left( C(t_2)/C(t_1)\right) /(t_2-t_1)\) and \(t_1<t_2\). Then, the infection doubling time is

$$\begin{aligned} \Delta T = (t_2-t_1)\frac{\log (2)}{\log \left( C(t_2)/C(t_1)\right) }. \end{aligned}$$

In this work, we use the average doubling time computed over consecutive days in the first month (four weeks) of the semester (see Fig. 2).

2.3 Global Sensitivity Analysis

In this work, we are interested in understanding how particular COVID-19 SEIR epidemic model factors \(\varvec{\theta }=(\theta _1,\theta _2,\dots ,\theta _k)\), Table 1, affect a model response Y (the cumulative number of infections, C(t), and the infection doubling time, \(\Delta T\)). We perform a variance-based global sensitivity analysis on these epidemic dynamics with the Sobol method (Saltelli et al. 2010, 2008; Archer et al. 1997), an approach to decompose the response variance by single and combined factor interactions

$$\begin{aligned} \text {Var}(Y) = \sum _{i} \text {Var}_{\theta _i}+\sum _i\sum _{j>i}\text {Var}_{\theta _{i,j}}+\cdots +\text {Var}_{\theta _{1,2,\dots ,k}}. \end{aligned}$$

Under this decomposition, the proportion of variance from a single model factor, the first-order index, can be written as

$$\begin{aligned} \mathbb {S}_i = \dfrac{\textrm{Var}_{\theta _i}(E_{\varvec{ \theta \sim i}}(Y\vert \theta _i))}{\textrm{Var}(Y)}, \end{aligned}$$
(3)

where E is expectation, \(\textrm{Var}\) is the variance, \(\theta _i\) is the \(i^\text {th}\) model factor, and \(\varvec{\theta _{\sim i}}\) indicates varying all factors except \(\theta _i\). In our work, a second quantity of interest is the total-order index, or total-effect index (Homma and Saltelli 1996), which in addition to the first-order index information, accounts for the additional contribution of a factor to the model variance from interaction effects with other model factors

$$\begin{aligned} \mathbb {S}_{T_i} = \frac{E_{\varvec{ \theta \sim i}}(\textrm{Var}_{\theta _i}(Y\vert \theta _{\sim i}))}{\textrm{Var}(Y)}. \end{aligned}$$
(4)

It follows that \(\mathbb {S}_{T_i}\ge \mathbb {S}_i\), and that a total-order index value of zero indicates that the model factor is non-influential. In our analysis, we estimate the first-order and total-order indices through numerical model solutions generated by sampling from the input factor space following Saltelli et al. (2010).

To estimate the first-order index, \(\mathbb {S}_i\) (Eq. 3), we use the estimator

$$\begin{aligned} \frac{1}{N}\sum _{j=1}^N f(\varvec{B})_j \left( f \left( \varvec{A_B^\text{(i) }} \right) _j-f(\varvec{A})_j \right) , \end{aligned}$$

for a particular model factor i, first presented in Saltelli et al. (2010). Given a model response \(Y = f(\cdot )\), and where \(\varvec{A}\) and \(\varvec{B}\) are \(N\times k\) matrices in which each row is a sampling of the k model factors (N samples per matrix). Matrices \(\varvec{A}\) and \(\varvec{B}\) are generated using the same model factor sampling method and are thus interchangeable but serve for bookkeeping when forming the k matrices denoted by \(\varvec{A_B^{(*)}}\). The matrix \(\varvec{A_B^{(i)}}\) is an \(N\times k\) matrix where column i comes from matrix \(\varvec{B}\) and all other \(k-1\) columns come from matrix \(\varvec{A}\).

The total-order index for a particular model factor i, \(\mathbb {S}_{T_i}\) (Eq. 4), is estimated with the estimator

$$\begin{aligned} \frac{1}{2N}\sum _{j=1}^N \left( f(\varvec{A})_j-f \left( \varvec{A_B^\text{(i) }} \right) _j \right) ^2. \end{aligned}$$

This estimator was first presented in Jansen (1999). For a more technical and detailed explanation of these estimators, one may consult Saltelli et al. (2010).

To adequately sample the multidimensional input factor space, when forming the matrices \(\varvec{A}\) and \(\varvec{B}\), we use Latin hypercube sampling implemented in MATLAB using the Statistics and Machine Learning Toolbox (MathWorks 2020). To compute the first-order and total-order indices, we compute \(N(k+2)=31,500\) model solutions (\(N(k+1)\) model solutions are suggested in Archer et al. (1997) to estimate Sobol indices) with \(N=1,500\) simulations for each of the \(k=19\) model factors. Each model factor is sampled from a uniform distribution over the range provided in Table 1. All confidence intervals (CIs) were formed by resampling the \(N(k+2)\) simulations 2000 times, with replacement, to compute the Sobol indices.

3 Data and Contacts

While the results we present here are specific to UC Merced, we note that any interested bubble-like communities could calculate their contacts as we outline below and use the model to analyze measures intended to mitigate COVID-19 transmission. For example, our model can be applied to skilled nursing facilities (where in-patients would be equivalent to ‘on-campus students,’ out-patients would represent ‘off campus students,’ doctors might correspond to ‘faculty/staff,’ and therapists/nurses might translate to ‘graduate students’ and ‘classes’ in this case would be face-to-face treatments). UC Merced is a public land-grant university set in a rural community. In these calculations, there were 8326 undergraduate students, of whom 2885 live on-campus and 5441 live off-campus. There were 721 graduate students and 424 faculty and staff. Lecturers and postdocs are considered part of the faculty/staff population. We expect our results would hold for similarly sized campuses with similar class structures. Although we are able to get exact data from the registrar, we provide a template that can be used to estimate the contacts between different subpopulations.

For our model to be accurate, it is necessary to be able to estimate the number of contacts each of the subpopulations (off-campus undergraduates, on-campus undergraduates, graduate students, and faculty/staff) has with one another. We assume that most of these contacts come from in-class instruction, living/dorm situations, contact with the outside community, and unscheduled social interactions. We explain each of these factors separately and denote \(\mathbb {C}^c\) for classroom contacts, \(\mathbb {C}^l\) for living/dorm contacts, \(\mathbb {C}^o\) for contact with the outside community, and \(\mathbb {C}^s\) for unscheduled social interactions. The sum of these four matrices provides the total contact matrix \(\mathbb {C}\). We use \(c_{i,j}^k\) to denote the (ij)-th entry of the submatrix \(\mathbb {C}^k\) with \(k\in \{c, l, o, s\}\). For simplicity, this information is displayed as a contact matrix in Table 2. We do not dynamically model the outside community and so the contact matrix is of size \(4\times 5\), since the interactions of campus population do not influence the outside community. The units of our contact matrix are contact-hours per day.

Table 2 The contact matrix \(\mathbb {C}\) represents the total contact-hours per day within and between subpopulations, broken down by classes (\(\mathbb {C}^c\)), living situations (\(\mathbb {C}^l\)), outside community engagement (\(\mathbb {C}^o\)), and unscheduled social interactions (\(\mathbb {C}^s\))

In the following subsections, we will discuss how the contact matrix is filled using in-class instruction, living situation, outside community interaction, and unscheduled social interaction information. Each of these sections corresponds to a different component of the contact matrix, as described in Table 2.

3.1 In-class Instruction \(\mathbb {C}^c\)

The majority of interactions in our model are derived from classroom instruction. We assume that lectures are comprised of 3 h per week, while discussion sections meet for one hour per week. When a class is listed as a laboratory, we assume that it meets for 2.5 h per week. Further, we assume that faculty and staff only interact with each other during faculty meetings. In particular, we assume a faculty meeting occurs once every other week for an hour. The ‘average’ department size was calculated using information from UC Merced School of Natural Sciences and was roughly 17.5 faculty per department.

Table 3 A schematic for calculating the contact due to teaching and class interactions at a universityy. Since classes only occur for the university community, this represents a \(4 \times 4\) submatrix of Table  2 In other words, let \(\mathbb {C}^{cs}\) denote the matrix calculated from the schematic shown here, then the full classroom contact matrix would be \(\mathbb {C}^c=[\mathbb {C}^{cs}, \textbf{0}]\) where \(\textbf{0}=(0;0;0;0)\)

Especially in large classrooms, counting every other individual as a contact just as likely to spread COVID-19 is not realistic. COVID-19 is spread through droplets (e.g., sneezing, coughing) for contacts within 6 feet, but it can also be spread through aerosols which remain suspended in the air and viable for much longer (Brlek et al. 2020; Hamner 2020; Kwon et al. 2020). As such, our formulation for the spread of COVID-19 in a classroom distinguishes between ‘close’ transmission (through droplets) and ‘far’ transmission (through aerosols). First, we assume that 8 contacts (the 8 closest people surrounding you) have a ‘normal’ chance of infecting you. Then, beyond your 8 neighbors, you have a reduced fraction (25%) of the ‘normal’ chance. (We also ran results with a 10% chance of infection, see Supplementary File 1, but the sensitivity results were qualitatively similar to the 25% chance.) As an example, assume an off-campus undergraduate student is in a 3-hour per week lecture course with 40 other undergraduate students, 10 of whom live on-campus and 30 of whom live off-campus, then their ‘contact-hours’ from this one class with other off-campus undergraduate students would be

$$\begin{aligned} \begin{array}{cc} \Bigg (\underbrace{6}_{\begin{array}{c} \text { close contacts} \\ \text {off-campus UG} \end{array}} + \underbrace{(30-6)}_{\begin{array}{c} \text { far contacts} \\ \text {off-campus UG} \end{array}}\times 0.25\Bigg )\times \underbrace{3/7,}_{\begin{array}{c} \text { hours spent in}\\ \text {this class per day} \end{array}} \end{array} \end{aligned}$$

and with on-campus undergraduate students would be

$$\begin{aligned} \begin{array}{cc} \Bigg (\underbrace{2}_{\begin{array}{c} \text { close contacts} \\ \text {on-campus UG} \end{array}} + \underbrace{(10-2)}_{\begin{array}{c} \text { far contacts} \\ \text {on-campus UG} \end{array}}\times 0.25\Bigg )\times \underbrace{3/7.}_{\begin{array}{c} \text { hours spent in}\\ \text {this class per day} \end{array}} \end{array} \end{aligned}$$

Note we assume that of the 8 contacts closest to the student, a proportional amount come from each sub-population. We compute these values for each of their courses and sum them together to get the total number for that student. We average over all the individuals in a subpopulation to obtain the corresponding classroom contact matrix entries.

We display the relevant information for calculating classroom contacts in Table 3. Classroom contacts only influence the campus subpopulations (not the outside community), and thus, the fifth column (not displayed in Table 3) consists of all zeros. We note that in Table 3, the classes graduate students teach as graduate assistants may be laboratories or discussion sections, the classes faculty and staff teach may be lectures, laboratories, or discussion sections. For our baseline contact matrices and more details, see Appendix A.

3.1.1 Network Analysis

Although we use the contact matrices to model the connectivity of the campus, we can also use network analysis since we have access to individual-level data. In particular, we can build a network graph of every individual on campus and their connectivity to each other via classes only. From this network, we can determine characteristics of those individuals that have the highest rates of contact. These network analyses are presented in Appendix B.

To build the weighted undirected graph, we consider each individual affiliated with the university as a separate node. Edges are formed between two nodes if those two individuals share a class (either as students or as student/instructor). For each additional class individuals share, the edge weight is increased by one. Of course, the graph changes depending on whether all classes meet in person, or whether there is a class capacity. Supplementary Fig. B1 (see Supplementary File 2) displays the resulting network for a full campus in which all classes meet in person. The network is colored by subpopulation (reddish purple corresponds to on-campus undergraduates, yellow represents off-campus undergraduates, sky blue for graduate students, and vermilion for faculty/staff), the size of each node represents the weighted degree, and the edge thickness reflects weights. It is apparent that under normal conditions most of the campus is connected with one another. However, as seen in Supplementary Fig. B4 (see Supplementary File 2), when classes with more than 25 students enrolled do not meet in person, much of the network becomes disconnected, i.e., having no contact with any other individuals.

With a resulting graph for the no intervention strategy that has over 9000 nodes and 1.6 million unique edges, it becomes necessary to use metrics to analyze exactly how the campus is connected. We report a histogram distribution of the weighted degree, a measurement calculating how many edges each node has. Histogram representations of the degree are displayed in Fig. 12. It is clear that capping in-person classes has a dramatic effect on the weighted degree, reducing maximal degree from 1,417 to 367. In fact, edges are reduced to \(\approx \) 0.13 million. Further, approximately 17% of campus individuals have 0 classroom contacts (no edges) when implementing a class cap strategy.

We can gather information about possible ‘super-spreaders’ by examining the individuals that have the highest degrees. Under the assumption of no interventions, it is clear the individuals with the highest degree are most often undergraduate students taking many introductory classes. For example, the individuals with the top 3 degree scores were all undergraduate students taking 4–5 lower level introductory classes (all living off-campus). When class caps are implemented, the structure of those with the highest interactions changes. Of the top 3 degree scores when imposing a class cap of 25, we have one lecturer teaching multiple labs and two graduate students that are taking a full load of classes and also teaching discussion sections and labs. This highlights how important it is to incorporate these subpopulations, who may serve as a vector for disease transmission, in our model.

3.2 Living Situations \(\mathbb {C}^{\varvec{l}}\)

In addition to classroom contacts, most members of the community will also interact with other individuals based on their living situation. We assume that the living contacts are based on living situation only. Therefore, there is no mixing between subpopulations in our model. For example, on-campus undergraduate students do not live with off-campus undergraduate students and vice versa. This means that in our contact matrix (Table 2), living contacts exist only on the diagonals.

We assume that on-campus undergraduate students have contacts for approximately 10 h per day with their direct roommate. To simulate the effects of encountering other individuals in their dorm, we assume that on-campus undergraduate students have 2 h of contact per day multiplied by the average number of beds per bathroom. For example, if there are roughly 4 students per bathroom, they would experience 8 h of additional contact per day.

The off-campus undergraduates are assumed to live with, on average, 3 off-campus undergraduate housemates. Since their living situation is likely larger than a dorm room, we assume less contact, at 4 h per day contact with each roommate. Graduate students have a similar situation, except that we assume they live with 1.5 other graduate housemates. Similarly, we assume 4 h of contact per day with each roommate for each graduate student. We assume that faculty and staff to do not live with other faculty and staff.

3.3 Contact with Outside Community \(\mathbb {C}^o\)

One of the most important aspects about a bubble-like community, from an infectious disease perspective, is that some key individuals have contact with the “outside community.” This is often overlooked in mathematical models, partly because the populations that interact with the outside world tend to be outnumbered by those contained fully in the bubble. However, these outside contacts cannot be ignored because they represent the potential for infection to infiltrate the “closed” community. These contacts are present in the fifth column of the contact matrix displayed in Table 2.

There are varying levels of contact with the surrounding community depending on which subpopulation a person is part of. We assume that there is little contact with the outside world if you live on-campus (1 h of contact per day). For off-campus undergraduate students and graduate students, the number of contacts is higher due to increased shopping, transportation, etc., at 5 h of contact per day. We assume that faculty and staff have the highest amount of contact with the outside community since many faculty live with families that are not affiliated with the university (15 h of contact per day). However, as these numbers are not directly produced from known data, we also include a parameter c, a multiplier in front of the community contact matrix, that we vary.

3.4 Unscheduled Social Interactions \(\mathbb {C}^s\)

One aspect of contact that has not yet been addressed is contact that occurs outside the classroom and living situation. In particular, we consider the effect of “unscheduled” social interaction in which members of the undergraduate student population meet for gatherings, unmasked, on a daily basis. Examples of these daily social interactions might include eating dinner with friends in the dining hall or forming an in-person study group for a course. In our contact matrix (Table 2), these are included in the undergraduate student populations (both on- and off campus).

We incorporate unscheduled social interactions in our model in two ways: first the daily weekday interactions described above and secondly an increase in social interaction during the weekend. To calculate the daily social interactions, we scale the classroom contacts by 30%. Our assumption is that the amount of contacts each student has socially will be similar to those experienced in the classroom, especially as some of our “unscheduled” social interactions include course study groups. The scaling value is chosen as studies show that the contact-hours in ‘other/social’ is roughly 30% of the ‘school’ contacts (Glass and Glass 2008; Leung et al. 2017). We assume that as class caps are implemented, unscheduled social contacts similarly decrease. Since many of these unscheduled social contacts are based on study groups arising from courses, unscheduled social contacts decrease as class caps are implemented. Thus, the social contact matrix \(\mathbb {C}^s\) is always 30% of the classroom contact matrix \(\mathbb {C}^c\). There is large uncertainty in number of social contacts, and so we multiply this contact matrix by the parameter p, which varies from 50 to 150% as shown in Table 1 as part of the sensitivity analysis. For the increased social interaction over the weekend, we multiply these daily social contacts by the parameter \(\omega \), to simulate going to a larger gathering. This parameter and larger number of contacts is active from 5 pm Friday until 5 pm Sunday.

4 Results: Global Sensitivity Analysis

We are interested in the sensitivity of the cumulative number of infections at a point in time since the semester began, \(C(t_i)\), and the sensitivity of the infection doubling time, \(\Delta T\), to the model parameters and initial conditions. We refer to initial conditions and model parameters as “model factors” in this text. We performed a global sensitivity analysis following Saltelli et al. (2010) and estimated the first-order sensitivity index and the total-order effect through numerical model solutions generated by sampling from the input factor space, assuming that all parameters were uniformly distributed in their given range listed in Table 1.

4.1 Variance in Cumulative Infections and Infection Doubling Time

We first study the variance of the following model metrics: the time-varying cumulative infections, the doubling time, and number of the cumulative infections at the end of the fifteen-week term. As mentioned above, we do this by employing a global sensitivity analysis approach where we vary parameters independently and uniformly over their ranges. In our analysis, we consider the quantifiable effect of three classes of model factors: infection parameters (\(\beta \), \(\sigma \), \(\phi \), \(\kappa ^R\), \(\kappa ^V\)), contact parameters (M, c, \(\omega \), p, m, \(\alpha \)), and initial conditions (\(I_{x,0}^{y}\) where \(y \in \{ s,a\}\) and \(x\in \{ d,u,g,f\}\)) (see Table 1 and Appendix A for model factor ranges and details).

Fig. 3
figure 3

(Color figure online) Distribution of cumulative infections. Figures show the distribution of cumulative infections over the span of a semester (15 weeks) and we allow the contact, infection parameters, and initial number of infectious individuals to vary (see Table 1). The mean and coefficient of variation for the doubling time and total cumulative infections are reported in each subplot. The reduction in expected cumulative infections (RECI) from case (a) with no vaccination and no class caps is presented in each subsequent case

Figure 3 presents the behavior of the cumulative number of infections as a consequence of class caps (each column represents a different class cap: none, 100 students and 50 students) and vaccination status of the campus at the beginning of the term (each row). For simplicity, we assume that all undergraduates have the same vaccination fraction \(v_{u,0}=v_{d,0}\) and that graduate students and faculty have the same vaccination fraction \(v_{f,0} = v_{g,0}\), at the start of the semester. We consider the following three vaccination scenarios, Low: \(v_{u,0}=v_{d,0} = 0\%, v_{f,0} = v_{g,0} = 0\%\); Medium: \(v_{u,0}=v_{d,0} = 40\%, v_{f,0} = v_{g,0}=50\%\); and High: \(v_{u,0}=v_{d,0} = 80\%, v_{f,0} = v_{g,0}= 100\%\). The black line signifies the mean cumulative infections over time, while the pink and blue shadings show one and two standard deviations from the mean, respectively. The mean and coefficient of variation for the doubling time and total cumulative infections are reported in each subplot.

Without vaccination or class caps (Fig. 3a), we expect about 1500 infections by the end of the semester (the baseline number of cumulative infections). By implementing a class cap of 50 students, the expected cumulative infections drop to about 200, an 86.8% reduction in cumulative infections (RECI) from the baseline (Fig. 3c). On the other hand, increasing vaccination rates to 100% for faculty and graduate students and 80% for undergraduate students reduces the expected cumulative infections to about 112 (92.6% RECI) by the end of the semester (Fig. 3g). With a class cap of 50 and high vaccination rates, we expect about 82 cumulative infections (94.6% RECI) by the end of the semester (Fig. 3i).

Fig. 4
figure 4

(Color figure online) Expected infection doubling times and cumulative infections by class capacity and percent of vaccinated undergraduates. The expected cumulative number of infection (cumulative infections) by the end of the semester and the expected doubling time computed during the first four weeks of the semester. The error bars are a 95% confidence interval

In addition to examining the cumulative number of infections, we can also investigate how vaccination and class caps impact the doubling times of infections among the campus populations. Figure 4 displays the infection doubling time in weeks (top row) and cumulative infections at the end of the semester (bottom row) as a function of vaccinated undergraduate students, \(v_{u,0} = v_{d,0}\). The columns correspond to faculty and graduate student vaccination rates 0% (left), 50% (middle), and 100% (right). Within each panel, the effects of incorporating class caps are portrayed with no class cap (yellow/thick line), 100-person cap (pink/medium line), and 50-person cap (teal/thin line). As we increase undergraduate student vaccination rates, we observe the lengthening of infection doubling time and we note that infection doubling times often exceed the length of the semester (Fig. 4a–c). We also see that increasing faculty and graduate student vaccination rates is not as effective in reducing the number of cumulative infections as vaccinating the larger undergraduate student population. This can be seen in Fig. 4d–f in the similarity between the three panels where the faculty vaccination is increased by 50% from left to right. The curves representing the cumulative infections, as a function of undergraduate vaccination rate, indicate a minimal change in the number of cumulative infections as faculty and graduate student vaccination rate is increased to 100%. This is likely because the faculty and graduate students only make up about 12% of the campus population.

It is apparent that the NPI of having large enrollment classes be remote drastically reduces the number of cumulative infections by the end of the semester, especially when the campus population is not significantly vaccinated (Fig. 4d–f). In the case when none of the population was vaccinated, by capping classes at 50 students, there is an 86.8% reduction in the expected cumulative infections by the end of the semester (Fig. 3c). Increasing the percentage of the vaccinated population also has a large effect in reducing the cumulative infections by the end of the semester. In the case when there is no class cap, having 80% of undergraduates and 100% of faculty, staff, and graduate students vaccinated resulted in a 92.6% reduction in the expected cumulative infections by the end of the semester (Fig. 3g). We observe from Fig. 4, that under low vaccination, class caps have a higher effect at reducing the number of cumulative infections. However, under high vaccination, implementing class caps has less of an effect in reducing the expected number of cumulative infections by the end of the semester.

4.2 Sobol Analysis of the Variance in Cumulative Infections and Infection Doubling Time

The first-order Sobol sensitivity index \(\mathbb {S}_i\) measures the direct effect that each model factor \(\theta _i\) \(\in (\beta \), \(\sigma \), \(\phi \), \(\kappa ^R\), \(\kappa ^V\), M, c, \(\omega \), p, m, \(\alpha \), \(I_{u,0}^s\), \(I_{u,0}^a\), \(I_{d,0}^s\), \(I_{d,0}^a\), \(I_{g,0}^s\), \(I_{g,0}^a\), \(I_{f,0}^s\), \(I_{f,0}^a)\) has on the variance of the model output. The total-order Sobol sensitivity index \(\mathbb {S}_{T_i}\) measures the total effect (direct and through interactions with other model factors) each model factor, \(\theta _i\), has on the variance of the model output. As mentioned above, we consider the sensitivity of two model outputs: the cumulative number of infections over time and the doubling time of the infection. We categorize the parameters for the sensitivity analysis into three groups: infection parameters (\(\beta \), \(\sigma \), \(\phi \), \(\kappa ^R\), \(\kappa ^V\)), contact parameters (M, c, \(\omega \), p, m, \(\alpha \)), and initial conditions (\(I_{u,0}^s\), \(I_{u,0}^a\), \(I_{d,0}^s\), \(I_{d,0}^a\), \(I_{g,0}^s\), \(I_{g,0}^a\), \(I_{f,0}^s\), \(I_{f,0}^a\)). In each sensitivity analysis figure, we show 9 subfigures illustrating three class cap scenarios—no class cap (the first column), class cap with 100 students (the second column), and class cap with 50 students (the third column)—along with three vaccination scenarios—(1) \(0\%\) vaccination (the top row); (2) \(50\%\) of faculty, \(50\%\) graduate students, and \(40\%\) of undergraduate students vaccinated (the middle row); and (3) \(100\%\) of faculty, \(100\%\) graduate students, and \(80\%\) of undergraduate students vaccinated (the bottom row).

We note that the initial conditions do not contribute to the sensitivity of either of our metrics in a fashion that is dependent on the vaccination or NPIs. As such, we include those figures in Appendix C.

4.2.1 Sensitivity of Infection Doubling Time

We now turn to examine the global sensitivity analysis with respect to infection doubling time \((\Delta T)\). Figure 5 displays the global sensitivity analysis of infection doubling time for infection and contact model parameters, while Fig. 13 presents the global sensitivity analysis of infection doubling time with respect to initial conditions.

Fig. 5
figure 5

(Color figure online) Global sensitivity analysis of infection and contact parameters on infection doubling time. First-order (blue) and total-order (red) Sobol indices are shown as well as the standard errors. The mean and coefficient of variation for the doubling time and total cumulative infections are reported in each subplot. Each column represents three class cap scenarios: none, 100 student, and 50 student caps. Each row represents one of three vaccination scenarios at the start of the semester. First row: \(0\%\) vaccination; second row: \(50\%\) of faculty, \(50\%\) graduate students, and \(40\%\) of undergraduate students vaccinated; third row: \(100\%\) of faculty, \(100\%\) graduate students, and \(80\%\) of undergraduate students vaccinated

When examining Fig. 5, we can see that across these vaccination scenarios and class caps, infection transmission rate \(\beta \) is the most significant parameter having the greatest individual effect on the variance in infection doubling time (\(\Delta T\)). Across these scenarios, transmission rate, \(\beta \), captures 40% to 60% of the variance in first-order effects, while considering interaction effects with other model factors (the total-order effect), \(\beta \) captures 40% to 70% of the variance in infection doubling time. That is, \(\beta \) has the largest effect on infection doubling time \(\Delta T\) among all model factors. We also see that both greater vaccination rates and lower class caps lead to a decrease in the effect of transmission rate on infection doubling time. The strong sensitivity of infection doubling time on transmission rate implies that anything that can be done to lower the transmission rate, such as wearing masks or improving HVAC systems, can have a large impact on infection dynamics at the beginning of the semester.

Other parameters to which infection doubling time is sensitive to include contact parameters: infections from the outside community, M, community contact multiplier, c, weekend contact multiplier, \(\omega \), social contact percentage, p, mask usage, m, and the probability that a symptomatic individual will self-isolate, \(\alpha \). From Fig. 5a, under no class cap and no vaccination we observe an infection doubling time of \(\Delta T = 3.5\) weeks (the number of cumulative infections would double in 3.5 weeks) and the sensitivity of infection doubling time is similarly affected by transmission from the community (through M and c) and transmission on campus (through p, \(\omega \), and m). The sensitivity from parameters of each type of contact, off-campus or on-campus contacts, capture about 10% of the variance in infection doubling time; however, as the class cap is reduced to 50 students, Fig. 5c, infection doubling time increases to \(\Delta T = 7.2\) weeks and the effect of mask, social, and weekend contact parameters is reduced. In this case, the effect of infections from the outside community, M, and community contact multiplier, c, becomes more influential. This indicates that when contacts are reduced on campus through NPIs that reduce the number of people that meet in one location, this reduces the effect from the transmission rate \(\beta \) and the social and weekend contacts, and infections from community and contact with the community have a greater effect on infection doubling time of cumulative on-campus infections. When vaccination is increased, as in Fig. 5g, under no class caps, with faculty and graduate student vaccination rates of 100%, and 80% of the undergraduate student population vaccinated, at the beginning of the term, we see an increase in infection doubling time to \(\Delta T = 16.1\) weeks (longer than our 15 week term) but no real qualitative change in the sensitivity of model factors when compared with reducing class cap. This minimal change in sensitivity to parameter values is because vaccination reduces the number of infected and susceptible individuals, but vaccination does not reduce contacts in the way that implementing class caps does by directly reducing person to person interactions. This analysis highlights the importance of including community interactions in models of bubble-like communities. Moreover, this shows that an effort to reduce transmission on campus will require attention to infection rates in the surrounding community.

Across scenarios in Fig. 5, the probability that a symptomatic individual will self-isolate, \(\alpha \), and the amount of time spent in the ‘exposed class,’ \(\sigma \), have nonzero total-order effect and no first-order effect. This indicates that \(\sigma \) and \(\alpha \) alone do not contribute significantly to the model variance, but through interaction with other parameters can have a significant impact on infection doubling time. Unfortunately, determining exactly which model factors \(\sigma \) and \(\alpha \) are interacting with, and to what extent, is not computationally tractable with this Sobol analysis. Another parameter that has a nonzero effect on doubling time is m or mask usage, where \(1-m\) is the reduction in infection transmission probability by wearing a mask. From Fig. 5, we observe that as we increase vaccination rates (moving down a column), we can see that the impact of wearing a mask, model factor m, stays rather consistent in its importance. This indicates that masks are still crucial for controlling the spread of disease even as vaccination rates increase. However, when decreasing the class caps, we can see that the sensitivity to m decreases. Thus, NPIs such as wearing masks and reducing contact by moving courses online are still important in controlling the number of infections on campus. The duration of naturally acquired immunity \(1/\kappa ^R\) and immunity through vaccination \(1/\kappa ^V\), do not affect infection doubling time since both types of immunity are assumed to last longer than three months (see Table 1) and infection doubling time is computed using cumulative infections observed within the first month.

Figure 13 displays the sensitivity of infection doubling time to the number of initially infectious individuals on campus. The distribution and number of initially infectious individuals at the start of the semester (undergraduate students, graduate students, and faculty and staff) do not play a large role in determining infection doubling time. This is intuitive, since infection doubling time measures the time it takes to double the number of cumulative infections and should be the same whether we start with one infectious individual or 100 infectious individuals (see Fig. 2).

Fig. 6
figure 6

(Color figure online) Global sensitivity analysis of infection and contact parameters on cumulative infections at the end of the term. First-order (blue) and total-order (red) Sobol indices are shown as well as the standard errors. The mean and coefficient of variation for the doubling time and total cumulative infections are reported in each subplot. Each column represents three class cap scenarios: none, 100 student, and 50 student caps. Each row represents one of three vaccination scenarios at the start of the semester. First row: \(0\%\) vaccination; second row: \(50\%\) of faculty, \(50\%\) graduate students, and \(40\%\) of undergraduate students vaccinated; third row: \(100\%\) of faculty, \(100\%\) graduate students, and \(80\%\) of undergraduate students vaccinated

4.2.2 Sensitivity of Cumulative Infections at End of Term

Here, we are interested in quantifying the effect of each parameter on the cumulative number of infections at the end of the term for each vaccination and class cap scenario. Figure 6 displays the global sensitivity analysis for the cumulative infections at the end of the semester with respect to infection and contact parameters. We find that as vaccination rates increase and class caps are implemented, the most sensitive parameters begin to change. Under the scenario with no vaccinations and no class caps (Fig. 6a), the most important parameters are the transmission rate, \(\beta \), weekend contact multiplier, \(\omega \), social contact percentage, p, mask usage, m, and the probability that a symptomatic individual decides to self-isolate, \(\alpha \). By incorporating different class caps, the landscape of sensitivity continues to change. The first-order effect of self-isolation, \(\alpha \), mask usage, m, social contact percentage, p, and weekend contact multiplier, \(\omega \) become less important, while the parameters governing the outside community start to play a larger role. The probability that Merced community individuals are infectious, M, and the community contact multiplier, c, start to influence the cumulative number of infections. This underscores the importance of educating local communities to help reduce the spread of infections in the community. As vaccination rates are increased, without class caps, we see a similar shift towards a decrease in first-order effects from mask usage, m, social contact percentage, p, and weekend contact multiplier \(\omega \); however, these parameters remain influential through interaction effects with other parameters. When both vaccination and class caps are utilized, the community parameters become more important than the transmission rate, see Fig. 6i.

We can also examine the sensitivity of cumulative infections to the initial number of infectious individuals on campus. Figure 14 presents this sensitivity analysis. It is apparent that under no interventions, the initial number of infectious individuals has little effect on the cumulative number of infected individuals at the end of the term. Similarly, when incorporating class caps, the initial conditions are insensitive. However, under high vaccination rates (80% undergraduate students and 100% faculty, staff, and graduate students), the cumulative number of infections is slightly sensitive to the initially infectious individuals on campus—with slightly higher sensitivity to undergraduate student infections, both symptomatic and asymptomatic. The sensitivity in those scenarios is present because the final number of cumulative infections observed at the end of the term is mainly from those infected at the start of the semester and undergraduate students make up the greatest proportion of the campus population (about 88% of the on-campus population consists of undergraduates). This tells us that under the best-case scenario with both high vaccination and class caps, it is important to minimize the number of infectious students at the start of the term.

4.2.3 Sensitivity of Time-Varying Cumulative Infections

In order to dissect the forces that contribute most strongly to the cumulative infections, we next consider their time-varying first and total-order Sobol indices. That is, we now look at the total contributions to the variance in cumulative infections over time. We separate the impact of infection parameters (Figs. 7, 9), contact parameters (Figs. 8, 10) and initial conditions (Fig. 15, 16). Because the time-varying first and total order indices for initial conditions show that their effect is limited to the start of the term, and that their effect is greatest when the initially infected made up the majority of cumulative infections observed in one semester, we include those figures in Appendix C.

Fig. 7
figure 7

(Color figure online) Time-varying total-order effect of infection parameters on cumulative infections. The total-order Sobol indices are shown as well as the standard errors. The mean and coefficient of variation for the doubling time and total cumulative infections are reported in each subplot. Each column represents three class cap scenarios: none, 100 student, and 50 student caps. Each row represents one of three vaccination scenarios at the start of the semester. First row: \(0\%\) vaccination; second row: \(50\%\) of faculty, \(50\%\) graduate students, and \(40\%\) of undergraduate students vaccinated; third row: \(100\%\) of faculty, \(100\%\) graduate students, and \(80\%\) of undergraduate students vaccinated

First, we consider the time varying impact of infection parameters on the cumulative infections. As shown in Fig. 7, 9, regardless of the vaccination status or class cap scenario, the strongest contributor at all times comes from \(\beta \). This is consistent with the results from Fig. 5, where we look at the sensitivity to the infection doubling time. We now see that this effect increases during the early part of the academic term. In the case of no vaccination and no class caps, Fig. 7a, the amount of time spent in the ‘exposed state,’ \(1/\sigma \), becomes influential early in the semester and reaches its strongest contribution to cumulative infections during the fourth week of the semester, the effect from \(\sigma \) then decreases. Both increasing vaccination and implementing class caps have the effect of delaying and dampening the contribution from \(\sigma \). The probability of becoming asymptomatic, \(\phi \), at most explains about 5% of the variance in the case with no vaccination and a class cap of 100 students. Also, the duration of immunity from both vaccination, \(1/\kappa ^V\), and natural immunity, \(1/\kappa ^R\), does not affect the variance in cumulative infections, including through interactions with other parameters. This is expected since we are only modeling a 15-week term, and both natural and immunity through vaccination lasts at least 12 weeks (see Table 1).

Fig. 8
figure 8

(Color figure online) Time-varying total-order effect of contact parameters on cumulative infections. Each column represents three class cap scenarios: none, 100 student, and 50 student caps. Each row represents one of three vaccination scenarios at the start of the semester. First row: \(0\%\) vaccination; second row: \(50\%\) of faculty, \(50\%\) graduate students, and \(40\%\) of undergraduate students vaccinated; third row: \(100\%\) of faculty, \(100\%\) graduate students, and \(80\%\) of undergraduate students vaccinated

Next, we present the role of contact parameters in Fig. 8, 10. In Fig. 8, we note that importance of contact parameters depends not only on vaccination and class cap scenario, but also on time. In the case of no vaccination and no class cap, Fig. 8a, we see that both the weekend contact multiplier, \(\omega \), and social contact percentage, p, are initially the most influential among all the contact parameters, but their effect decreases after the sixth week. We observe that as vaccination is increased without decreasing the class cap, Figs. 8d, g, the effect on cumulative infections from these parameters is delayed towards later in the semester. However, as the person-to-person contacts on campus are reduced through class caps as in Fig. 8b, c, we see the sensitivity to \(\omega \) and p replaced by the sensitivity to M, the percentage of infectious people in the community and c, the community contact multiplier. This is similar to the results from the sensitivity analysis in Figs. 5, 6. In the case with no vaccination and no class cap, Fig. 8a, we can see how \(\alpha \), the probability that a symptomatic person will self-isolate, becomes more influential over time but reaches a fixed effect of about 10% of the variance. Without class caps, as vaccination is increased to medium and high vaccination, Fig. 8d, g, respectively, we see that the effect due to \(\alpha \) is delayed to later in the semester but still reaches 10% of the variance by the end of the semester. However, under low vaccination and a class cap of 100 students (Fig. 8b), \(\alpha \) becomes more significant, explaining more than 15% of the variance by the end of the semester. When comparing Fig. 8, the total-order effect, and Fig. 10, the first-order index, in terms of \(\alpha \) we can see that most of the contribution to the model variance seen in Fig. 8 comes from the interaction effects with other model parameters. This means that in bubble-like communities with low vaccination rates and no class caps policy, making sure that students self-isolate if they are symptomatic plays a large role in affecting the cumulative number of infections observed over the course of the semester.

Fig. 9
figure 9

(Color figure online) Time-varying first-order effect of infection parameters on cumulative infections. Each column represents three class cap scenarios: none, 100 student, and 50 student caps. Each row represents one of three vaccination scenarios at the start of the semester. First row: \(0\%\) vaccination; second row: \(50\%\) of faculty, \(50\%\) graduate students, and \(40\%\) of undergraduate students vaccinated; third row: \(100\%\) of faculty, \(100\%\) graduate students, and \(80\%\) of undergraduate students vaccinated

Fig. 10
figure 10

(Color figure online) Time-varying first-order effect of contact parameters on cumulative infections. Each column represents three class cap scenarios: none, 100 student, and 50 student caps. Each row represents one of three vaccination scenarios at the start of the semester. First row: \(0\%\) vaccination; second row: \(50\%\) of faculty, \(50\%\) graduate students, and \(40\%\) of undergraduate students vaccinated; third row: \(100\%\) of faculty, \(100\%\) graduate students, and \(80\%\) of undergraduate students vaccinated

From Figs. 15, 16, we see that the infectious people at the beginning of the term only affect the cumulative infections within the first few weeks of the term. While we do observe different levels of sensitivity to these initial conditions, these effects are separated by population size. Among those that are infectious at the start of the term, those that are symptomatic have a slightly larger effect than those that are asymptomatic even though asymptomatic people are 50% less infectious than symptomatic people in our model. The largest effect on cumulative infections, from up to 1% of the population that is initially infectious, comes from undergraduate students that live off-campus (5441 students), followed by students living in the dorms (2885 students) and the least influential infections, consistently explaining less than 1% of the variance, are both graduate student (721 graduate students) and faculty/staff initially infected (424 faculty and staff). In all cases, we see that both increased vaccination and implementing lower class caps extend how far into the term these model factors remain significant. For example, in Figs. 15i, 16i, we see that the effect from the initially infectious undergraduate populations persists until the end of the term. That is because high vaccination and smaller class caps leads to fewer cumulative infections by the end of the term, and the initial infections make up a large number of the cumulative infections observed over the term.

For nonlinear models, it is not generally true that the first and total order Sobol indices will be equal. This only happens in cases where the model factors only have a linear contribution to the model variance as in the case of the initial conditions (Fig. 15, 16). However, in this case we note that the importance of model factors (i.e., the rank of their contribution to the variance) is similar for first-order and total order time varying Sobol indices (compare Figs. 7, 8, 9 and 10). However, there are a few intriguing differences to point out. Most notably, the first-order index for \(\beta \) (Fig. 9a) may exhibit an internal peak rather than simply increase to saturation. Our interpretation again would involve a change in the dynamics from early in the semester with a large susceptible population to later in the semester with less susceptible people that can become infected. For the case with no class cap and no vaccination, the transmission rate \(\beta \) is the parameter that has the most significant direct effect on the output variance (Figs. 9a, 10a). This indicates that any interventions that can lower the transmission rate, \(\beta \), would significantly impede the transmission. As shown in Fig. 9a, we observe that the fraction of variance due to \(\beta \) begins to decrease after week 3 and begins increasing again after week 6. From Fig. 7a of the time varying total-order effect, the saturation in \(\beta \) is likely due to the increasing interaction effects with the contact parameters (Fig.8), and with the other model factors that have a nonzero effect.

5 Discussion and Conclusion

In this manuscript, we introduced an ODE-based SEIR model with multiple subpopulations (undergraduate students living on-campus, undergraduate students living off-campus, graduate students, and faculty and staff) on a campus that interacts with the outside community. We discussed how to use registrar information to estimate contact-hours due to classes. We examined the effects of NPIs such as social distancing in the form of transitioning large classes to an online format, mask usage, and considered differing vaccination rates with waning immunity for both natural immunity and immunity from vaccination among the different campus subpopulations. We acknowledge that there are other methods of incorporating heterogeneous interactions, such as agent-based modeling and stochastic modeling (Bahl et al. 2021; Panovska-Griffiths et al. 2022). Our proposed model preserves some aspect of heterogeneity while remaining open to traditional methods of analysis and being computationally inexpensive. This allowed us to evaluate different intervention strategies and obtain rapid results to assist the administration in making decisions regarding the mitigation of COVID-19 on the university campus.

We examined the sensitivity of two epidemic characteristics, the infection doubling time and the cumulative number of infections over the course of the semester, to 19 model factors that fall into three categories: infection parameters, contact parameters, and initial conditions, under varying levels of vaccination rates and class caps. Through our analysis, we found that implementing class caps (or limiting the number of people that can meet in an enclosed space) and greater vaccination rates, led to changes in the rank importance of the model factors in their contribution to the model variance. We found that, in scenarios with increased vaccination rates without class caps, the social contact multiplier, weekend contact multiplier, and probability to self-isolation when symptomatic are the most significant model factors contributing to the infection doubling time and the number of cumulative infections observed over the course of one 15-week academic term (one semester). However, when decreasing person-to-person contacts through smaller class caps, under low vaccination, the contacts with infectious people from the outside community play a predominant role in affecting the number of cumulative infections for the campus population and directly affects the infection doubling time. This work highlights that even when the campus implements and regulates NPIs, or can require vaccination, the outside community still have a large effect on the infections observed on campus. Therefore, it will be necessary for universities to work with their surrounding communities to help limit the spread of COVID-19.

Many of the parameters, such as the contact matrices, were directly informed from campus data. Thus, it may be that these results may not hold for other campuses that differ from our own (e.g., no large classes, campuses that are more integrated with the surrounding community, or campuses in urban settings). However, our flexible framework allows other universities to alter parameters and calculate contact matrices from their registrar data in order to make decisions about intervention strategies and to perform sensitivity analyses for their own campuses.

The parameters that were not directly calculated from data remain a source of uncertainty in the model. Although many of the parameter values were taken from literature values, the accuracy of those estimates is unknown. Therefore, model validation with real data remains a future direction for exploration. We also plan to replace the way that we model infections originating from the surrounding community, which is currently a constant, with dynamic data from county dashboard reports. We can also study the impact of a surge in COVID-19 infections in the off-campus community on the campus population by introducing a pulse in the community parameter.

We note there are several limitations with the current study. In our model, it is assumed that vaccination makes the recipient 100% immune with the possibility of breakthrough cases through waning immunity. A logical future direction would include booster shots and the introduction of novel variants (e.g., delta, omicron). This will be especially important when modeling an academic year as opposed to one 15-week academic term. Another modeling option would be to include vaccinated subpopulations with decreased transmission rates and the possibility of shorter convalescence times. As differing variants have become dominant, the parameters would need to be updated to reflect appropriate transmission rates and, perhaps, a shortened infectious period of time. We also assumed 100% compliance with interventions, meaning that every member of the campus community wears well-fitting masks at all times when they are indoors. In our work, we attempted to capture the varying degrees of mask usage and effectiveness by varying the mask efficiency from having no effect to reducing transmission by up to 60%.

Overall, our model exhibits results consistent with current public health messaging. Our model suggests that, despite increasing vaccination rates, it remains important to continue to socially distance, to wear masks, and to implement class caps to reduce the transmission of COVID-19. It is likely that these results will not hold as more COVID-19 variants emerge; however, our modeling framework can be readily adapted to account for such novel infections.