The coronavirus was first identified in Wuhan, a city of over 11 million people in China's Hubei province, in December of 2019. COVID-19, the disease caused by the coronavirus, grew by several thousand per day in China between late January and the peak of their epidemic in early February. The number of infections appearing each day has decreased there significantly. This was mainly due to stringent containment efforts.
However, the outbreak is now a global pandemic. Large outbreaks in South Korea, Iran, Italy, Spain, and the United States have all sparked interest in coronavirus research. As of April 19th 2020, there have been 2,432,092 cases and 166,794 deaths confirmed worldwide. Also, the World Health Organization has officially named COVID-19 a global pandemic. The spread of the virus has been rapid, with 208 countries now having reported cases.
Lawmakers and public health experts such as National Institute of Allergy and Infectious Diseases director, Dr. Anthony Fauci, have been encouraging "social distancing" as a measure to flatten the curve. This implies that we want to push down the peak number of infected people and spread it out over time. Many people understand how flattening the curve impacts the healthcare system from an empirical or qualitative perspective, see Goebel [8] . However, there is little quantitative research in this area and relatively little research on how much pushing is needed to flatten the curve.
In this paper, we model the potential demand for hospital resources like beds, ventilators, medical staff as an infinite server queue with non-stationary arrival rates. We use as our simple examples of unimodal arrival rate functions, the scaled Gaussian distribution density and the scaled gamma density functions. We then derive explicit expressions for the peak mean number in the queueing system as well as the time of this peak.
Our paper has two main goals. The first one is showing the global community how operations research, represented here by queueing theory, helps us understand the impact of COVID-19 on the healthcare system. The other goal is showing that by leveraging insights from queues with time varying rates leads to simple descriptions of the hospital patient dynamics for the COVID-19 pandemic. More specifically, we show using a scaled Gaussian arrival rate and a scaled Gamma arrival rate, how to calculate the time of the peak infection rate, the peak mean number of infected people, and the time between the peak hospitalization rate and the time of the peak number of patients. We use an infinite server queueing model to demonstrate the impact of social distancing and flattening the curves of arrival rates and mean counts of on hospital capacity.
Operationally, we can learn from the situations we have already observed in China, Europe, and New York City to avoid similar outcomes in other places? From the data, we see that "every hotspot has its own curve" see Hilk et al. [11] . One important advantage of our work is that it can be used to study each individual curve for individual regions around the world. The geographic and temporal clustering of outbreaks have the potential to create a health care system collapse. Many states in the United States have imposed "stay at home" orders to contain the spread of the disease. However, since many people are asymptomatic, we can show through the analysis of our queueing examples that containment measures implemented in recent days may take weeks to have an observable effect in the data.
Although hospital bed capacity is an important concern for health care officials, it is not the most crucial bottleneck for providing critical care when patients need it. After observing many deaths in Italy, many experts and government officials are concerned about the number of available ventilators. Patients who become critically ill need individual ventilators but they have become a scarce resource. Equally important resources are medical specialists and staff members that can supervise the patients in their beds and operate the ventilators properly. Fundamentally, this paper strives to help healthcare managers understand what supply is needed for all these lifesaving resources at the peak of the ultimate curve for their aggregate demand.
The contributions of this work can be summarized as follows:
• Given a non-stationary arrival rate, we derive an expression for the mean number of infected individuals using an infinite server queue.
• Using the infinite server queue model, we derive the exact time of the peak load and the value of the peak load.
• Using our peak load calculations, we determine the value of "flattening" that needs to be done in order to make the peak load infection value below a pre-specified capacity level.
• We study the impact of the duration of the virus on the peak load dynamics of hospital resources and how it affects our flattening policies.
• We characterize the nonlinear relationship between flattening the curve and the time lag between the peak rate of newly admitted patients and the peak demand for hospital resources. We also show how the lag can be observed in real data.
We model the number of patients as a M t /G/∞ queueing model. There are two important reasons to begin with the M t /G/∞ queueing model even though it has an infinite number of servers and this does not seem realistic from a practical perspective. The first reason is that the M t /G/∞ queue is very tractable since the distribution is known in closed form for any non-stationary arrival rate. The second reason is that the M t /G/∞ queue is the best type of queue one can hope for where everyone is immediately served and no one ever waits for service. From the perspective of the COVID-19 epidemic, this means that anyone who wants a bed gets a bed immediately or who ever needs a ventilator gets a ventilator immediately. In this regard, the M t /G/∞ infinite server queue is an upper bound for queueing models with a finite number of servers and without abandonment since it represents the dynamics if the manager had access to an infinite amount of resources and is not resource constrained. Moreover, the infinite server queue has a history of being used to staff finite server systems, see for example Jennings et al. [12] , Feldman et al. [7] , Liu and Whitt [15] , Massey and Pender [16] . Thus, the infinite server queue serves as the first step to understanding more complicated systems.
In this section, we state the closed form formulas for the M t /G/∞ queueing model, which exploits the results of Eick et al. [5, 6] for the time varying infinite server queue. In the paper of Eick et al. [5] , they use the properties of the Poisson arrival process and use Poisson random measure arguments to show that the M t /G/∞ queue Q ∞ (t), has a Poisson distribution with time varying mean q ∞ (t). As observed in Eick et al. [5, 6] , q ∞ (t) has the following integral representation
where λ(u) is the time varying arrival rate and S represents a service time with distribution G, G = 1 − G(t) = P(S > t), and S e is a random variable with distribution that follows the stationary excess of residual-lifetime cdf G e , defined by
(2.5)
Remark. As a point of emphasis, we would like to remind readers that the word curve in this paper maps to two different things: the arrival rate function and the mean offered load. The general public might call these the infection rate and the number of patients in the hospital.
Moreover, there is a causal relationship between the two curves, i.e. you need to flatten the first one in order to flatten the second one.
In this section, we describe the dynamics of the M t /G/∞ queue whose arrival rate is driven by a Gaussian distribution function. Thus, the arrival rate function λ(t) is given by the following expression
where τ is mean of the Gaussian distribution, σ is the standard deviation of the Gaussian distribution, λ * is the number of total expected number of infected individuals over the lifetime of the epidemic, and ϕ (x) is the probability density function of the standard Gaussian distribution i.e.
One of the important things to observe is that the standard Gaussian pdf is bounded from above i.e.
(2.8) Figure 1 : Plot of Gaussian arrival rate functions.
Using the Gaussian pdf as our arrival rate gives us the following integral relationships.
In Figure 1 , we plot two Gaussian arrival rate functions. The blue curve is a Gaussian(τ = 10, σ = 4) and the red curve (flattened curve) is a Gaussian(τ = 20, σ = 4). By increasing the standard deviation by a factor of 2, we have reduced the peak value of the arrival rate by a factor of two. This follows from the mode of Gaussian being inversely proportional to the standard deviation σ. Now that we have a full understanding of the Gaussian arrival rate function that we will use for patient arrivals, we can leverage insights from non-stationary queues to understand the dynamics of the total number infected and how long they stay infected.
Theorem 2.1. The M t /G/∞ queueing model with a Gaussian distribution arrival rate has a Poisson(q M/G/∞ (t)) distribution where the mean q M/G/∞ (t) is the solution to the following ordinary differential equation
and the solution is given by
Moreover, for any value of t, we have that
Proof. We actually start with the solution of the queue length given in Equation 2.13. The solution is easily given by the formula from Eick et al. [5] since we have that
(2.14)
To find the differential equation, one simply takes the derivative of the solution with respect to the time parameter t. The reason that we get the final expression in Equation 2.13 is that the derivative of the Gaussian distribution yields the first Hermite polynomial times the Gaussian pdf. Finally, the bound on the standard Gaussian density function yields the bounds on the queue length. This completes the proof. Table 1 : Comparison between regular and flattened curves peak values and times (λ * = 100).
In Figure 2 , we plot the same two Gaussian arrival functions given in Figure 1 . However, to the right of each arrival curve, we also plot the subsequent queue length when an exponential distribution is assumed for the service distribution. The parameters for each curve and the service distributions are given in Table 1 . Table 1 and Figure 2 yield many observations.
The first observation is that we see that both arrival rate curves intersect with the queue length at the peak queue length. This is a nice property and can be explained by first noticing that E[S] = 1/µ = 1 in our example. Thus, at the peak queue length, we know the time derivative of the queue length is equal to zero i.e.
(2. 15) This implies that at the time of the peak queue length, t * we have
A second observation is that for this value of the service distribution, E[S] = 1, the peak arrival rate is larger than the peak queue length. We find that for the flattened arrival curve, the value of the peak queue length is closer to the value of the peak arrival rate than in the non-flattened case. This implies that even though one might flatten the curve and reduce the peak arrival rate by one half, it does not mean that the peak queue length are also reduced by one half. We observe that the queue peak queue length is 91% of the peak arrival rate in the non-flattened curve, while it is 97% of the flattened curve.
In addition to finding the solution of the queue length, we can compute an expression for the time of the peak queue length. We know that the time of the peak arrival rate occurs at the mode τ and our goal is to understand the lag effect i.e. the length of time the peak queue length occurs after the peak arrival rate occurs. Theorem 2.2 suggests that the time of the peak queue length occurs after the peak arrival time τ . This is generally well known in the queueing literature. However, this is an important fact for the COVID-19 outbreak as it helps determine exactly from the peak number of infections, when the peak load will be seen in the healthcare system.
Theorem 2.2. The time of the peak queue length is the solution to the following fixed point equation
Proof. The proof is found in the Appendix. Table 2 : Comparison between regular and flattened curves peak values and times (λ * = 100). Table 3 : Comparison between regular and flattened curves peak values and times (λ * = 100). .
In Figures 3 -4 , we plot our Gaussian arrival curves with their respective queue lengths. In Figure 3 , we set the mean service time to be E[S] = 2 and in Figure 4 , we set the mean service time to be E[S] = 10. On the left of the both Figures 3 -4 , we plot the actual arrival rate and the actual queue length; however ,on the right of Figures 3 -4, we plot the actual queue length, but an arrival rate that is scaled by the mean service time E[S]. We observe that the in the unscaled plots on the left, the arrival rate does not intersect with the queue length at the peak. However, in the scaled plots on the right, we observe that the arrival rate intersects with the queue length curve at the queue length peak. This verifies Equation 2.16 again. Finally, we observe in Table 1 that by reducing the arrival rate peak by one half reduces the queue length peak by one half when E[S] = 1. However, we observe in Tables 2 -3 that as we increase E[S], we see a smaller reduction in the peak queue length. Certainly we expect the infectious period to be longer than one day for COVID-19 so our theoretical results imply that we may need to do more flattening of the curve as patients spend more time infected.
Theorem 2.3. The time of peak queue length is bounded by the following upper and lower bounds
Proof. First we observe that from Theorem 2.2 that t * solves the following fixed point equation
Now to prove the upper bound on the peak time, one observes that the numerator is bounded by the standard Gaussian bound given in Equation 2.39. For the lower bound, one just recognizes that both expectations E S e · ϕ t * −τ −Se σ and E ϕ t * −τ −Se σ are non-negative. To make an equivalence with the last equality, we use the equality between the expected stationary excess random variable and the second moment of the service distribution i.e.
. This completes the proof. Theorem 2.3 provides upper and lower bounds on the time of the peak queue length. Numerically, one can use a bisection algorithm in order to explicitly determine the time of peak infection. Once the time of the peak infection is known, then one can use Theorem 2.1 to determine the queue length value at that peak time as well as the value of the peak infection. This provides estimates on the number of ventilators or beds needed at the peak time of infection.
Theorem 2.4. Let q M/G/∞ (t) be a M t /G/∞ queueing model with a Gaussian distribution arrival rate. If C is the maximum capacity queue length for all values of t, then the value of σ * that flattens the peak queue length below C is given by
(2.20)
Proof. We know the queue length is bounded by the constant
Now we let the constant equal the capacity and solve for the standard deviation to get
(2.23)
In this section, we specify the service distribution to be a constant. Although the constant assumption for the service time is quite unrealistic, this will simplify many of the complicated equations and yield simple insights. Our first result shows the queue length can be described by an ode and the solution is simply an integral of the arrival rate evaluated at the current time and the current time minus the constant service rate. Finally, we also show that the lag between the peak arrival rate and peak load is equal to half of the service time.
Corollary 2.5. The M t /D/∞ queueing model with a Gaussian distribution arrival rate is the solution to the following ordinary differential equation
and the solution is given by
Finally, the time of the peak is equal to
(2.28)
Proof. This follows from standard analysis and the symmetry of the Gaussian distribution.
When the service time distribution is exponential, the mean queue length, q ∞ (t), solves the autonomous ordinary differential equation
(2.29)
The subsequent ordinary differential equation is linear and has a closed form solution given by the following expression q ∞ (t) = q 0 e −µt + e −µt t 0 λ(s)e µs ds.
(2.30) Corollary 2.6. The M t /M/∞ queue with arrival rate given by a Gaussian distribution with parameters has the following closed form expression for the mean transient queue length
Theorem 2.7. When the arrival rate is a Gaussian function and the service time is given by an exponential distribution, then the lag in the time between the peak arrival rate and the time of the peak queue length is given by the following expression
Proof. The proof is given in the Appendix.
Thus, for the exponential distribution, we have an explicit formula for the lag in terms of the ψ(·) function. This function is related to conditional expectations of Gaussian random variables and is related to the inverse Mills ratio i.e
In the context of operations research, the ψ(·) function has been analyzed in resource sharing settings in Hampshire et al. [9, 10] , as well as, in the context of using risk measures for staffing multi-server queueing systems in Pender [21] .
In this section, we describe the dynamics of the M t /G/∞ queue whose arrival rate is driven by a Gamma distribution function. Thus, the arrival rate function λ(t) is given by the following expression
where α is shape parameter of the Gamma distribution, β is the rate parameter Gamma distribution, λ * is the parameter. Unlike the standard Gaussian, the Gamma distribution is not guaranteed to be bounded unless the parameter α ≥ 1. Thus, for the remainder of this paper, we assume that α ≥ 1. Under this condition, we have that
Using the Gamma pdf as our arrival rate gives us the following integral relationships
One reason for using the Gamma distribution is that it is not symmetric like the Gaussian distribution. From some of the data that we see in practice, the arrival rate distribution does not appear to be symmetric and may have a heavy tail. For example in Figure ? ?, we see that the tail of the data is heavier than the Gaussian distribution tail. This fact of the data encourages us to explore the Gamma arrival rate function which is asymmetric.
In Figure 5 , we plot two gamma arrival rate functions. The blue curve is a Gamma(α = 5, β = .5) and the red curve (flattened curve) is a Gamma(α = 10, σ = .5). By increasing the variance by a factor of 2, we have reduced the peak value of the arrival rate by a factor of two. This follows from the mode of gamma being inversely proportional to the standard deviation σ. Now that we have a full understanding of the gamma arrival rate function that we will use for patient arrivals, we can leverage insights from non-stationary queues to understand the dynamics of the total number infected and how long they stay infected.
Theorem 2.8. The M t /G/∞ queueing model with a Gamma(α, β) distribution arrival rate has a Poisson(q M/G/∞ (t)) distribution where the mean q M/G/∞ (t) is the solution to the following ordinary differential equation
and the solution is given by
Moreover, for any value of t and α ≥ 1, we have that
Proof. We actually start with the solution. The solution is easily given by the the formula from Eick et al. [5] . To find the differential equation, one simply takes the derivative of the solution with respect to the time parameter t. Finally, the bound on the standard Gaussian density function yields the bounds on the queue length.
Theorem 2.9. The time of the peak queue length is the solution to the following fixed point equation
.
Proof. In order to find the peak queue length we need to set the time derivative of the queue length to zero i.e.
• q (t * ) = 0.
(2.45) Now using the differential equation given in Equation 2 .41, we have
Isolating t * by itself on the left hand side we obtain
.
This completes the proof. What is nice about this result is that it decomposes into two parts. The first part is the mode of the arrival rate function, which in the case of the Gamma distribution occurs at α−1 β . The second part is the positive shift from the mode of the arrival rate function, which is
. This also implies that the peak queue length occurs after the peak arrival rate like in the Gaussian case. 
In this section, we use some of the data on COVID-19 infections and deaths, which was made available through Johns Hopkins University website Dong et al. [4] . In Figure 7 , we plot the number of confirmed infections per day in the countries of China, South Korea, Italy, and the United States. One immediately can tell these countries apart since the confirmed cases start earlier in China, then South Korea, then Italy, and then finally the United States. We also observe that as of April 16th, it appears that the United States has hit its peak and is on its way down. On each of the plots, we also overlay a Gaussian and gamma distribution function to get a sense what distribution parameters best fit the data. For China, the best Gaussian fit was a Gaussian(λ * = 314101, τ = 23, σ = 8.28) and the best gamma fit was a Gamma(λ * = 369218, α = 6.75, β = .25). For South Korea, the best Gaussian fit was a Gaussian(λ * = 34914, τ = 42, σ = 16.37) and the best gamma fit was a Gamma(λ * = 31220, α = 9.4, β = .2). For Italy, the best Gaussian fit was a Gaussian(λ * = 440393, τ = 60, σ = 26.79) and the best gamma fit was a Gamma(λ * = 340126, α = 9.57, β = .15). For the United States, the best Gaussian fit was a Gaussian(λ * = 1295233, τ = 80, σ = 14.72) and the best gamma fit was a Gamma(λ * = 1768739, α = 17, β = .2). In Figure 8 , we plot the arrivals and deaths per day from January 22nd to April 16th normalized by their maximum value (arrival = black dots, deaths = red dots). For China, South Korea, Italy, and the United States, the maximum confirmed infections was (15, 136) , (851), (6557), and (35,098) respectively. Moreover, for China, South Korea, Italy, and the United States, the maximum confirmed deaths was (252), (11), (919), and (4591) respectively. We observe that for most all countries except Italy, the largest number of deaths lags the largest number of confirmed infections. This is the lag effect that we observe in our queueing models. This lag effect also inspires our construction of a new way to view the lag effect.
In Figure 9 , we plot the cumulative distribution function of the total number of arrivals and deaths for the countries of China, South Korea, Italy, United States, Japan, and Germany. This plot matches the quantiles of the arrivals to the quantiles to that of the deaths. In this plot we observe for 5 out of the 6 countries that the cdf of the death count lags behind the arrival cdf. The only country that does not follow the lag effect is Japan and this observation is consistent with Japan not testing its citizens enough Montanari [19] . Moreover, we also observe that the lag effect is most pronounced in South Korea, which has done a great job of testing its citizens and keeping them alive Kim [14] . If we rank the countries according to their lag effects, it is from best to worst (South Korea, Germany, China, United States, Italy, Japan). This method of viewing the arrival and death cdfs allows one to understand the impact of testing and the quality of the healthcare system in that country. These observations are also confirmed in Table 6 where we plot the lag according to the quantile levels in the cdf plots given in Figure 9 . Country q = .1 q = .2 q = .3 q = .4 q = .5 q = .6 q = .7 q = .8 q = .9  China  2  4  5  6  6  6  9  11  14  Germany  9  9  8  9  8  7  7  6  4  Italy  4  3  3  3  3  3  2  2  2  Japan  -4  -7  -5  -2  -2  2  2  1  1  South Korea  5  9  13  18  20  21  22  18  12  United States  5  5  5  5  4  3  3  3  2   Table 6 : Lags between arrival and death cdfs by quantile.
The analysis of the queueing model provides both qualitative insights and prescriptive guidance on how to flatten the demand for hospital resources.
First, we show that the peak number of hospitalized patients is inversely related to flattening the rate of arriving patients (Theorem 2.1). This insight provides a mechanism through the standard deviation parameter σ to control the peak number of hospitalized people. However, the effectiveness of social distancing on flattening the curve, i.e. increasing σ, remains an open question. As a result, in order to avoid overcrowding flattening should be inversely proportional to the amount of hospital capacity. Additionally, we find that flattening the curve has a smaller impact on the peak number of patients when the service times are long. This means the longer patients stay in the hospital, the more society must flatten the curve in order to achieve the same capacity level. In the three examples with an exponential service distribution and a Gaussian arrival rate, we flattened the peak arrival rate by 50%, however, we see that the peak queue length was reduced by 47%, 42%, 23% respectively for E[S] = {1,2,10}. For the gamma arrival rate setting, we see that a 33% reduction in the peak arrival rate, we obtain a 31% and 16% reduction respectively for E[S] = {1,10}.
A key takeaway from the queueing analysis is that the time of the peak arrival rate of infected patients does not coincide with the time of the peak of number of hospitalized people (or deaths). There is a time lag. We characterize this lag exactly for the simplified queue model (Theorem 2.2). First, we show that the peak number of infected lags behind the peak number of new infected individuals. Surprisingly, there is a nonlinear relationship between flattening the curve and the time lag between the peak rate of newly admitted patients and the peak demand for hospital resources (Theorem 2.2). Moreover, we show in Theorem 6.2 that the lag, , is bounded by the mean service time for an exponential distribution i.e.
This analysis also implies that as we flatten the curve i.e. increase σ, the lag converges to the mean service time 1 µ .
In this paper, we present a simple infinite server queueing model of the number of infected people with Covid-19. This analysis can be used to understand how many people will become infected, need a hospital bed, or need a ventilator. We show how the dynamics of the number infected explicitly depends on the arrival rate function. We explore two arrival rate functions. The first is symmetric and is given by a Gaussian distribution function and the second is is asymmetric and is given by a gamma distribution function. We explicitly calculate the time of the peak queue length, the peak queue length, and its dependence on the duration of infection. Despite the insights derived in this paper, this is not the whole story and many questions still remain. In what follows, we outline some of the important questions that still remain and need to be studied. By leveraging operations research tools and techniques, hospitals have prided themselves on the efficiency that comes with just-in-time supply management and minimizing empty beds, thus saving on hospital costs. However as a result, they are often not equipped for an epidemic surge such as this one. Despite our efforts to understand the impact of COVID-19 on the healthcare system, there is plenty of additional research that needs to be done in order to understand emerging complexities. One area of importance is to understand the effects of surges in arrivals into the healthcare system. Some recent research has studied queueing models with self-exciting point processes that model random surges in demand, see for example Daw and Pender [1, 2, 3] , Pender [21] , Niyirora and Pender [20] , Massey and Pender [18] . These types of models would generalize our model in this paper, but, they are more complicated and would disguise the easy to understand insights from simpler models.
In this work we have not explicitly analyzed the impact of testing and lags in testing for the coronavirus. It has been shown through empirical analyses that testing delays can have a huge impact on the reported numbers of positive cases and health outcomes, see for example Ryan-Mosley [23] . We believe that more research in understanding how testing can affect the spread of the disease is important. It is clear from our cdf plots in Figure 9 that testing is important for predicting deaths that are to follow.
Finally, we are observing in some cities that it is important to use demographic information for understanding the impact of COVID-19 on particular communities or regions of the country. In the city of Chicago, it has been observed that even though African Americans make up roughly 30% of Chicago's population, African Americans make up 52% of cases of COVID-19 and 68% of deaths, see for example Reyes et al. [22] . Also in NYC, it has been reported that COVID-19 is impacting low-income communities much harder than those that are well off Wilson [24] . Recent work by Kahara and Pender [13] combines queueing theory with demographic information to get a deeper understanding of blood donation dynamics and it is clear that similar demographic analyses are also needed to understand the full impact of COVID-19 on marginalized populations.
Proof. First we observe that from Theorem 2.2 that t * solves the following fixed point equation
Thus, if we define X to be a unit exponential random variable, then we have
µσ · e −µ · e µ 2 2 /2 · ( − µσ 2 ) · Φ σ − µσ + σ · ϕ σ − µσ µσ · e −µ · e µ 2 2 /2 · Φ σ − µσ (6.62) Now it remains to prove the following equalities E ϕ − X /µ σ = µσ · e −µ · e µ 2 2 /2 · Φ σ − µσ (6.66)
For the first equality we define G to be a Gaussian random variable and we will exploit