Controlling the spread of COVID-19 on college campuses

1.   Introduction 

As most readers know, the COVID-19 pandemic began in China in December 2019, then slowly spread around the world. Among its impacts on the US are changes to the education system. After Spring Break 2020, which occurred in mid-March, many colleges sent students home and switched to remote learning. During the summer, the big question became what to do in the fall. Colleges across the US have moved classes online, restricted the number of students on campus, and have modified semester schedules to run from mid-August to Thanksgiving, which occurs on the fourth Thursday in November, with no breaks. To further control the spread of the virus on campus, they have made dining halls take-out only and reconfigured public spaces to encourage social distancing. In addition, there will be daily symptom reporting, increased testing, as well as quarantine and contact tracing for infected individuals.

In this paper we will consider two methods to suppress the spread of the virus. In Section 2, we ask the question: How much would spread be reduced if all students in dorms have single rooms? This is a significant issue for Duke because undergraduates are required to live on campus three of their four years, so at any time more than 80% live in dorms. By using a generalization of the household model ^[1], we show that the reduction can be substantial. Independent of the research in this paper Duke reduced the on campus student population so that every student could have a single room.

In Section 3, we turn our attention to classrooms and ask a two part question inspired by Gressman and Peck ^[2], who analyzed a complex stochastic agent based model to determine the impact of having large classes be online-only. The first part is: How does the distribution of class sizes affect virus spread? To approach this question, we considered two scenarios for a college that has $1000$ students, each taking $3$ classes. In Scenario 1, all classes have size $30$. In Scenario $2$, 1/4 of the classes have size 60, and 3/4 have size 20. By calculating the maximum eigenvalue of the matrix that governs the spread in the second case, see ^[3], we show that the basic reproduction number $R_0$ is 48% larger in scenario 2 due to increased spread in larger classes.

Knowing that larger classes promote virus spread, it is natural to move classes above a certain threshold online. Thus, the second question we ask in Section 3 is: Where should the threshold be for moving classes online? To address this question, we consider a scenario where 2405 students take three classes each and there is one class of each size from $10$ through $120$. This is not a very realistic situation but it does allow us to investigate a wide range of cutoffs. We find that the threshold for moving classes online to reduce the basic reproduction number $R_0 = 1$ and hence stop the spread of the epidemic depends on the basic infection probability p in the model. For the values we investigate, the threshold can be as large as 70 or as small as 30.

Here we have examined two aspects of campus life in isolation to have models that are mathematically tractable. ^[4,5] have also developed models that integrate dorms and classes with dining, gyms, and socialization. The household model is appropriate for our situation since all the students are roughly the same age. If we were instead modelling the spread through Durham, we would need an age structured mode with a matrix of rates that govern the interactions between the different age classes. See ^[6,7,8]. For spread through a state or country, one must consider spatial aspects ^[9,10,11]. For more about epidemics in general see ^[12]. For discussions of epidemics on random graphs see ^[13,14,15].

2.   Household model and college dorms 

The question that motivates this section is: To what extent would the spread of COVID-19 on the Duke campus be reduced if all students had single rooms? To answer this question, we will use results that have been developed for the household model. See ^[1] and references therein. The household model improves on homogeneous mixing models by partitioning the population into multiple households with the probability of infecting an individual being larger for individuals in the same household.

We will review the necessary theory on the household model in Section 2.1. In Section 2.2, we examine the case where all students have single rooms by viewing dorms as households and applying the household model. In Section 2.3, where we examine double rooms, we still view dorms as households, but generalize our approach on computing the dynamics within a household.

2.1.   Household model 

Let the population of size $N = mn$ be split into $m$ households of size $n$. For simplicity, we consider the case of fixed infection time. Let $p_H$ be the probability an individual $i$ infects individual $j$ who is another member of the same house. Let $p_G$ be the probability of infecting individual $j$ who is not in the same house as individual $i$. $H$ is used to denote household, and $G$ is used to denote global.

Step 1: Connect each pair of individuals in the same house with probability $p_H$.

Let $\pi_k$ be the probability an individual belongs to a cluster of size $k$. Let $\mu = \sum_k k \pi_k$ be the mean cluster size, and let $G_H(z) = \sum_{k = 1}^\infty \pi_k z^k$ be the generating function of the cluster size. In our analyses, we suppose that the epidemic is subcritical in each dorm. Results in ^[1] can be used to analyze the case in which the epidemic in a single dorm is supercritical. There are $> 40$ dorms on the Duke campus, so in this case there would be no hope of controlling the spread on campus. Fortunately, the single dorm supercritical scenario is not consistent with estimates of $R_0$ for COVID-19.

Step 2: Make long distance connections between individuals in different households with probability $p_G$.

The number of long distance infections caused by one person is $Y_1\approx$ Poisson($Np_G$) if $N$ is large. The generating function of $Y_1$ is $G_{\nu}(z) = \exp(Np_G(z-1))$.

From steps 1 and 2, we see that

To compute the generating function, we recall that if $X, Y\ge 0$ are independent with generating functions $\phi_X$ and $\phi_Y$, then $X+Y$ has generating function $\phi_X\phi_Y$. This implies that the generating function of the number of individuals infected by a cluster of size $k$ is $\exp(Np_G(z-1))^k$ and hence that the generating function for the number of individuals infected by the cluster of a randomly chosen individual is

A well-known result for branching processes implies that the probability $\rho = 1-\zeta$ that the process dies out is the fixed point $< 1$ of the generating unction of the offspring distribution:

Our epidemic is essentailly percolation on the graph of connections in the household model, so the solution $\zeta$ of Eq (2.1) is an excellent approximation of the probability of a large epidemic when $n$ and $N$ are large. The size of a large epidemic is the same as the size of the giant component in percolation, which is the same in our homogeneous graph to the probability of a large epidemic.

When we apply the household model to our dorm models by viewing a dorm as a household, $G_H$ will vary between models since the choice of single or double rooms results in a slightly altered local spread within a dorm. On the other hand, since the global spread of the infection is unaffected by room size, $G_\nu$ will be the same for the two dorm scenarios.

2.2.   Single rooms 

In adapting the household model to the spread of infection in a college dorm, $p_H$ will instead be written as $p_D$, since we are concerned with a dorm rather than a household. Similarly, $G_H$ will instead be written as $G_D^1$, with the superscript 1 indicating single rooms.

We first compute $R_0$. Consider $m$ dorms with $n$ single rooms, where population size $N = mn$. If $n$ is large, e.g., 120 is the size of a typical freshman dorm at Duke, then the total number of individuals infected in one dorm as a result of one infected is the total progeny of a branching process. The branching process will have the offspring distribution Poisson($\lambda_1$), where $\lambda_1 = np_D$. Here $n$ should be $n-1$ but we will ignore these small imprecisions in what follows. We assume $np_D < 1$ so that the mean of the total progeny is

Each individual infected in the dorm infects an average of $Np_G$ individuals outside the dorm so the basic reproduction number is

Now, we solve for $G^1_D$. Breaking things down according to the number of individuals directly infected by the first individual, the generating function of the size of the epidemic within the dorm caused by one infected is

where $G^1_b(z) = \exp(np_D(z-1))$, the generating function of the spread of the infection within the dorm. One method for solving Eq (2.4) is through iteration. We start with $h_0(z) = z$, which corresponds to 1 individual in generation 0, and then iterate

$h_{\ell}(z)$ is the generating function of the number of individuals infected in the first $\ell$ generations of the process, so $h_{\ell}(z) \downarrow h(z)$, a solution of the recursion.

In the case of single rooms, an explicit solution for $G^1_D$ can be found by taking the derivative of equation Eq (2.4) to get the differential equation

with the boundary condition $G^1_D(1) = 1$. It is well known that the solution to the differential equation is

where $W$ is the Lambert $W$ function which is defined by $W^{-1}(z) = ze^z$.

By Eq (2.1), the probability of a large epidemic is the largest solution of

This can be further simplified. Plugging in $G^1_D$ yields

Multiplying through by $-\lambda_1$, then applying $W^{-1}$ to both sides gives

Recalling that $\lambda_1 = np_D$ and simplifying, we obtain

which is the survival probability for a branching process with a Poisson($Np_G +np_D$) distribution.

2.3.   Double rooms 

Suppose now that we have $m_1$ dorms with $n_1$ double rooms, so that $N = 2m_1n_1$. For the sake of comparison, we let the number of people in a dorm be the same, regardless of whether they have only single rooms or only double rooms, i.e., $n = 2n_1$. Let $p_L$ be the probability of infecting your roommate, and $p_D$ be the probability of infecting an individual within your dorm. To find the analogue of $G^1_b$ in this scenario, we think of a generation as having a first step where the ancestor splits into two with probability $p_L$, stays 1 with probability $1-p_L$ and then a second step in which long range connections are formed within the dorm. The offspring distribution is now a mixture of two Poissons:

The mean number of infections in one generation, $\lambda_2$, will be the sum of the means of the two Poisson distributions so $\lambda_2 = (1+p_L)2n_1p_D$.

Imitating calculations in Section 2.2, if $\lambda_2 < 1$, the mean of the total progeny of the branching process is

and the basic reproduction number is

To compute $G^2_D$, the analogue of $G^1_D$ for the double room scenario, we note that the generating function of the number of infections in one generation is

so $G^2_D$ satisfies

$G^2_D$ can then be solved using the iteration method in Section 2.2.

Having computed $G^2_D$ we can compute the probability of a large epidemic as in the previous case by numerically solving

A comparison of the epidemic probabilities $\zeta$ for double-room and single-room dorms is given in Figure 3. Note that in the double-room situation when $2n_1p_D = 0.7$ or 0.9, there is a positive probability of a large epidemic within a single dorm, so the curves are positive when $Np_G = 0$. When $2n_1p_D = 0.3$ or 0.5, $p_G$ needs to exceed a positive threshold for there to be an epidemic in a university with double-room dorms, but those thresholds are higher when there are only single-room dorm.

3.   Spread of infection in college classes 

This section addresses two questions. The first asks how different class size distributions affect pandemic spread. The second asks whether taking classes above a certain size online can prevent outbreaks.

To address these questions, we consider some simplified models in which each student takes three classes. The situation can be described by a graph in which there are $n$ students and $m$ classes with sizes $c_i \; 1\le i \le m$, and there is an edge from each student to the three classes they are enrolled in. Each classroom is assumed to be homogeneously mixing. To formulate the dynamics, we say that the three class members that correspond to one student are their alter egos.

An individual that is infected in class $i$ at time, adds his two alter egos (his presence in other classes) to the infected population at time $1/2$ and then infections at time $1$ are produced. None of the infecteds from time $1/2$ are present at time $1$.

3.1.   Equal class sizes reduce spread 

To show that class sizes impact spread, we compare two scenarios for a college with $1000$ students that each take three classes.

Scenario 1. 100 classes of size 30. An individual infected in class $1$ at time 0 will produce a Poisson($29p$) number of new infected individuals in class $1$ at time 1. His two alter egos are equally likely to be in classes 2–100, and each alter ego will produce Poisson($29p$) infecteds in whichever class they are in. So if $m_{i, j}$ is the mean number of infecteds in class $j$ at time 1 from a single infected person in class $i$, then

From this we see that an infected at time 0 produces an average of $3\cdot 29p$ infecteds at time 1 so

Scenario 2. 25 classes of size 60 ($1 \le i \le 25$) and 75 of size 20 ($26 \le i \le 100$). Let $m(i, j)$ be the expected number of infections in class $j$ caused by one infected student in $i$. $c_i$ and $c_j$ represent that sizes of class i and j respectively. In this scenario, the probability that an individual in class $i$ has an alter ego in a different class $j$ is $2c_j/(3000-c_i)$, where $c_i$ is the size of class $i$. The $2$ comes from $2$ alter egos. Thus, we have

Since the eigenvector $x$ associated with the maximum eigenvalue will have $x_i = a$ for $1\le i \le 25$ and $x_i = b$ for $26\le i \le 100$, the eigenvalue problem for the $100 \times 100$ matrix can be reduced to the eigenvalue problem for the $2\times 2$ matrix with entries

Solving a quadratic shows that this matrix has eigenvalues $129.380p$ and $25.288p$. Thus,

This value is 48% larger than the $R_0$ for scenario 1, and suggests that larger classes increase virus transmission.

3.2.   Stochastic model for scenarios 1 and 2 

For a more detailed look at the two scenarios, we turn to a stochastic model. As in the previous section the model takes place in discrete time. Each time step is roughly one class meeting, so that two time steps are roughly equivalent to a week in real time. An infected person in a class infects each other individual with probability $p-0.01$ and has a 0.5 probability of entering quarantine. Using Eq (3.1) we see that $R_0 = 2(0.87) = 1.74$ in scenario 1. Eq (3.2) implies that $R_0 = 2(1.2938) = 2.5876$. $R_0 = 2.5$ is a commonly used value for an uncontrolled COVID-19 epidemic. Figures 4–6 confirm that Scenario 2 is worse than Scenario 1.

The next figure compares the equal class size and 60–20 scenarios. In both scenarios, the outbreak generally grows out of control within a few weeks and infects the majority of students. In Scenario 2, the outbreaks tend to develop faster, have a higher peak, and vary less in size and dynamics. Scenario 2 also produces a higher probability a large outbreak will occur from one individual, and larger final epidemic sizes. Thus, it is clear that the presence of larger classes, even if average class size remains the same, increases the risks for students on campus.

3.3.   Moving classes online 

In this section, we consider the effects of moving classes above a given size online. To do so, we consider a university with 111 classes, one each of sizes ranging from 10, 11... 120. The total class size is $10+11+...+120 = 7215$ and there are $7215/3 = 2405$ students. To begin, we first consider when all classes are in person.

Following the calculations for scenario 2 in Section 3.1, we see that the average number for people infected in classroom $j$ by an infected individual in classroom $i$ is

Similar to the previous scenarios in this section, we find the largest eigenvalue of the above matrix to determine $R_0$. Although we were unable to find the eigenvalue analytically, we used Python to determine that

When $p = 0.01$, $R_0 = 2.515$.

We now look at the situation where the university moves all of its classes of size $> k$. Viruses can't spread through online classes, we have

The other values of $m(i, j)$ remain unchanged. The change in $R_0$ depending on $k$ and various $p$ can be seen in Figure 7, which makes clear that establishing a cutoff can sharply reduce virus spread through classes. Although Gressman and Peck suggested a cutoff of $30$ based on their stochastic simulations for spread in universities ^[2], our results in Figure 7 suggest that the optimal cutoff depends largely on $p$. For sufficiently small $p$, even cutoffs as high as $70$ can be effective.

Acknowledgments 

The research for this paper was supported by the summer 2020 program DOmath, sponsored by the Duke Mathematics Department with partial funding from the Rhodes Information Initiative at Duke and Duke's Trinity Research Enhancement fund. Durrett and one of the undergraduate students were supported by NSF grant DMS 1614838 from the Math Biology Program.

Conflict of interest 

The authors have no conflicts of interest to disclose.