key: cord-0626467-im6njz64
authors: Awasthi, Achal; Minin, Vladimir; Huang, Jenny; Chow, Daniel; Xu, Jason
title: Fitting a Stochastic Model of Intensive Care Occupancy to Noisy Hospitalization Time Series
date: 2022-03-01
journal: nan
DOI: nan
sha: 5c8fe8693a34fd2316c4db8268ffd95fee48238e
doc_id: 626467
cord_uid: im6njz64

Intensive care occupancy is an important indicator of health care stress that has been used to guide policy decisions during the COVID-19 pandemic. To reliably forecast this occupancy, the rate at which patients are admitted to hospitals and intensive care units (ICUs) and the rate at which patients leave ICUs are crucial. Since individual-level hospital data are rarely available to modelers in each geographic locality of interest, it is important to be able to estimate these rates from publicly available daily numbers of hospital and ICU beds occupied. We develop such an estimation approach based on an immigration-death process that models fluctuations of ICU occupancy. Our flexible framework allows for immigration and death rates to depend on covariates, such as hospital bed occupancy and daily SARS-CoV-2 test positivity rate, which may drive changes in hospital ICU operations. We demonstrate via simulation studies that the proposed method performs well on noisy time series data, and apply our statistical framework to hospitalization data from the University of California, Irvine (UCI) Health and Orange County, California. By introducing a likelihood-based framework where immigration and death rates can vary with covariates, we find through rigorous model selection that hospitalization and positivity rates are crucial covariates for modeling ICU stay dynamics, and validate our per-patient ICU stay estimates using anonymized patient-level UCI hospital data.

In March 2020, the World Health Organization (WHO) declared the spread of a 2019 novel coronavirus (SARS-CoV-2) a global pandemic. Currently, we have seen over 350 million cases 1 , putting a tremendous strain on limited hospital capacities in countries around the world. Available hospital beds have proven to be one of the most valuable resources through the course of the pandemic, and in order to appropriately allocate resources, public health officials must rely on model forecasts to predict occupancy during the current and future pandemics [1] . Here, we propose a likelihoodbased method for estimating the rates at which coronavirus disease 2019 (COVID-19) patients enter and leave hospital intensive care units (ICUs) based on publicly available hospitalization time series.

To quantify the uncertainties surrounding novel diseases such as COVID-19, constant policy changes, and changing rates of vaccine distribution, stochastic epidemic models are well-suited compared to their deterministic counterparts.

However, stochastic modeling of SARS-CoV-2 transmission dynamics has not been more widely adopted, in large part due to computational burdens associated with inference for model classes such as partially observed Markov processes [2] . COVID-19 observational data are often incomplete, marked with coarseness and systematic underreporting. This often gives rise to intractable marginal likelihood functions associated with partially observed data.

Although likelihood-based methods provide principled ways to estimate parameters and quantify uncertainty, fitting a stochastic model in such missing data settings requires marginalizing over the set of all possible events that are compatible with the observed data. When this latent space is enormous, fitting stochastic epidemic models becomes a highly non-trivial exercise due to the increased computational complexity. This has spurred recent work on numerical approaches and approximation methods toward likelihood-based inference [3] , [4] , [5] , [6] , [7] .

In this article, we overcome this intractability for a class of immigration-death processes representing ICU occupancy.

We focus our attention on the simpler problem of estimating the dynamics and mean duration of ICU stays instead of the full dynamics of SARS-CoV-2 transmission through a population. Our approach accounts for the missing data probabilistically in a principled way by directly utilizing the likelihood of the observed data. We present a methodological framework for maximum likelihood estimation for a class of continuous time immigration-death processes from population-level data. Our models fall within the class of continuous-time Markov chains, and we carefully balance their flexibility with tractability to allow for deriving expressions for the transition probabilities of the process.

These comprise the marginal likelihood function; their solutions together with numerical and optimization techniques enable classical procedures such as maximum likelihood estimation. As a result, we can perform principled inference of interpretable model parameters in a mechanistic model, integral to understanding the mean duration of ICU stay by COVID-19 patients according to the underlying dynamic of the pandemic. Specifically, we begin by assuming that the per-patient immigration rate, i.e., the rate at which patients are admitted to the ICU, is proportional to the daily number of hospital beds occupied, while the per-patient death rate, i.e., the rate at which patients are discharged from the ICU, is constant over time. As previous studies have shown that the length of 1 https://coronavirus.jhu.edu/map.html hospital ICU stay varies depending on resource availability, patient characteristics, and hospital-specific differences in clinical care strategies (Leclerc et al., 2021), we extend the methodology to allow the model parameters to depend on a number of covaiates. The proposed framework is flexible as a result, allowing immigration and death rates to vary with a number of salient variables in a log-linear fashion.

The rest of the paper is structured as follows. In section 2.1, we present a description of the data that motivated us and also describe an immigration-death process in detail. In section 2.2, we introduce the the parametrization of the immigration rates based on relevant covariates and present our immigration-death framework in detail. In section 2.3, we present derivations of the Kolmogorov forward equations, transition probabilities, and the likelihood function for discretely observed data from an immigration-death process. We also briefly describe a technique that we use to calculate the transition probabilities in a computationally efficient way. In section 2.4, we describe the optimization methods used to numerically calculate the MLE and uncertainty quantification for the MLE.

In section 3.1, we design a simulation study to test our immigration-death framework. We then apply our model to hospitalization data from the UC Irvine hospital system as well as Orange County in sections 3.2 and 3.3 respectively, exploring the division of data into three phases following natural waves of the pandemic. The concluding section 4 discusses limitations and interpretations of the model as well as directions towards future work.

We begin by describing the data that motivate the immigration-death model. In our application to the problem of ICU stay estimation, the immigration and death rates will be interpretable as the rates at which patients are admitted to and discharged from the ICU respectively. In order to provide enough realism in the model without overfitting, we test and select among models of varying complexity. These methods hinge on access to the observed data likelihood, which additionally yields a natural approach to uncertainty quantification of the parameter estimates.

We consider hospital-level data collected from March, 2020 to December 2020, consisting of daily counts of occupied hospital beds and ICU beds in the University of California, Irvine's hospital system (UCI Health) and the larger Orange County (OC) level. Note that this surveillance level data does not contain any information about the patients at an individual level. The OC data are publicly available 2 , consisting of daily counts H 0 , H 1 , . . . H T of the number of beds occupied in OC hospitals by COVID-19 patients; we similarly denote by I 0 , I 1 , . . . I T the number of beds occupied daily in an ICU by COVID-19 patients. Plots visualizing the raw data appear in Appendix C. For the rest of the paper, hospitalization data refers to both the number of hospital beds occupied and the number of ICU beds occupied in a hospital.

Our goal is to fit a stochastic immigration-death model of ICU capacity to these data, inferring the per-patient rates of ICU admission as well as the rate of discharge. An immigration-death process X(t) is a continuous-time Markov chain (CTMC) defined over a countable state space of the non-negative integers. The process can increase by one in the case of an immigration event, decrease by one in the event of a death, or remain at the same state at any instant in time t > 0. The process X(t) can be completely characterized by its immigration rate α(t) and its death rate iµ(t)-we model the latter as proportional to the current state of the process, making the death process ratelinear. If the process is in state X(t) = i at time t, then the process can be defined by its instantaneous rates λ ij (t) which parametrize the infinitesimal generator matrix Λ Λ Λ(t) = {λ ij (t)} characterizing the time inhomogeneous CTMC.

To understand the relation of instantaneous rates to the behavior of the process, denote the transition probabilities P ij (s) = P(X(t + s) = j|X(t) = i), defined as the probability of the process transitioning from state i at time t to state j at time t + s. Then for a small time interval ∆t, the probabilities of an immigration event, a death event, and no event occurring are given by

As we model the ICU capacity as the immigration-death process X(t), naturally immigration can be interpreted as the flow of patients into the ICU, while death is interpreted as a patient leaving the ICU when either they are discharged or they die. We therefore view the count data {I 0 , I 1 , . . . I T } as discrete observations of the continuous-time process,

i.e. I 1 = X(t 1 ) for a sequence of observation times t i at which the counts are available. The immigration rate α(t)

will depend on hospitalization as well as other factors, described in the next subsection, while the per-person discharge rate µ(t) = µ is assumed to be homogeneous over the entire observation period [0, T ]. The summary of the model is illustrated in Figure 1 appearing in Appendix C. Note that the overall rate of ICU discharge iµ fluctuates proportionally with the current ICU occupancy, a natural modeling assumption in our context.

As the model described above may lack realism in positing that the immigration and (per-patient) death rates are constant throughout the entire observation period, we add flexibility by allowing them to depend on potentially timevarying covariates. We parametrize the log immigration rate as a linear combination of salient covariates while assuming a constant death rate over the entire observation period. Equation 2 describes such an immigration rate where log α plays the role of an intercept term, θ i (t) represents an i th covariate of interest at time t, and β i is its corresponding coefficient:

This parametrization enables us to account for a set of covariates that may be relevant in accurately describing how the model rates vary. As a large inflow of patients to the ICU comes from patients already admitted in the hospital, the number of patients admitted to the hospital is the leading predictor affecting the immigration rate in the above parametrization. We hypothesize that the proportion of infected people in a particular region is also likely to influence the rate of patients entering the ICUs. Thus, an overwhelming number of patients already admitted to the hospital combined with a higher proportion of infected people in a particular region suggest that the patients will enter the ICU at a higher rate. In particular, we allow the per-patient immigration rate into the ICU to depend on factors such as the daily number of hospital beds occupied H O (t), and the daily cumulative test positivity rate P A (t) in the region of interest.

The daily cumulative test positivity rate P A (t) is defined as the ratio of the cumulative number of individuals who Though it is intuitive that the potential covariates we consider would impact the rates of the process, it is preferable to quantitatively select the features to include in a way that balances model flexibility and generalizability given limited data. Within a likelihood-based framework, information-theoretic criteria allow for rigorous model selection. For instance, we may select covariates while guarding against overfitting using the Akaike information criterion (AIC)

where K is the number of free parameters and l(Y; Θ) is the maximized log-likelihood of the given data. Among a set of candidate models, the one with lowest AIC is preferred. We explore several competing models that posit different log-linear functional forms for α(t), summarized in Table 1 , where α, µ, β HP , and β H are non-negative coefficients. Table 1 : Various models considered that differ by the way they parametrize (log) immigration rates

# Parameters 1 log α(t) = log α + β H log H O (t) + β HP P A (t) log H O (t) + β P P A (t) 5 2 log α(t) = log α + β H log H O (t) + β HP P A (t) log H O (t) 4 3 log α(t) = log α + β H log H O (t) + β P P A (t) 4 4 log α(t) = log α + β H log H O (t) 3 5 log α(t) = log α + log H O (t) 2

Toward enabling likelihood-based inference, this subsection derives the transition probabilities of a simple immigration-death process, describes how the transition probabilities are related to the marginal likelihood function, and describes a technique to obtain the transition probabilities in a computationally efficient way.

Recall the transition probability P ij (s) = P(X(t + s) = j|X(t) = i) denotes the probability of the immigration-death process transitioning from state i at time t to state j at time t + s. Evaluating the likelihood in the proposed model framework relies on access to these quantities as they appear as a product in the likelihood function of the observed data.

We outline a method that derives their solution by way of the probability generating function (PGF), and also describe an efficient numerical way to compute them. This will allow us to numerically optimize the likelihood function toward parameter estimation. Similarly, we then detail how to perform interval estimation and covariate-based model selection efficiently using these numerical techniques.

Let {X(1), X(2) . . . X(T )} denote the numbers of occupied ICU beds at observation times t 1 , t 2 , . . . t T . By the Markov property, the likelihood function of these discretely observed data is given by a product of the transition

where Θ denotes the set of parameters we wish to infer. Thus, it is necessary to calculate the transition probabilities as they form the backbone of the observed likelihood function. Next we derive the transition probabilities for a discretely observed immigration-death process with immigration rate α(t) and death rate µ(t).

We begin by deriving the Kolmogorov's equations [ [8] ] for a discretely observed immigration-death process:

where the third equality is due to Equation 1. After subtracting P ij (t) from both the sides, dividing by h, and sending h → 0, we obtain the Kolmogorov equations:

Next we consider the PGF, 

subject to the initial condition φ i (s, 0) = s i . A detailed derivation of the steps toward obtaining and solving Kendall's PDE is available in Appendix A. The solution of the above PDE is

and provide a point of departure to calculate the transition probabilities. Notice the transition probability P ij (t) = P(X(t) = j|X(0) = i) can formally be obtained by differentiating the PGF φ i (s, t) and normalizing by an appropriate constant,

However, in practice repeated numerical differentiation is a computationally intensive procedure and often becomes numerically unstable, especially for a large j. Instead, we use the technique described in ( [10] ) and make use of the Fast Fourier transform (FFT) to calculate the transition probabilities P ij (t) which occur as coefficients of s j in Equation (4).

We map the domain s ∈ [0, 1] onto the boundary of the complex unit circle, considering a change of variables s = e 2πiω , so that the new domain becomes an equally spaced set of points along the unit circle in the complex plane. This allows us to view the generating function defined in Equation (4) as a periodic function

where the coefficients of interest P ij are now interpretable as the j th Fourier coefficients of a Fourier series. This suggests they can be recovered by the Fourier inversion formula,

A Riemann sum approximation can be used to discretize the integral,

with a larger choice of N yielding a more precise approximation. This approach avoids the need to directly compute numerical higher-order partial derivatives of φ i (s, t) as described in Equation (6). Instead, the FFT provides a fast way to extract all of the coefficients P ij (t) simultaneously for 0 ≤ j ≤ N − 1. (3) that the likelihood function of the discretely observed data is comprised of these transition probabilities. This FFT method therefore enables its efficient computation for any given set of parameter values, which can then be embedded within any generic optimization scheme that requires calls to the likelihood L(Y; Θ) as the objective function.

With an efficient way to compute transition probabilities in hand, we can use Equation (3) to numerically evaluate the likelihood function L(Y; Θ) for any parameter configuration, and in turn use generic iterative optimization methods to maximize it [11] . We find that the Nelder-Mead method [12] delivers robust yet efficient solutions to maximizing the log-likelihood, by way of a simplex algorithm using only calls to the objective function. That is, the method does not require gradient or other higher-order information [13] . Because we cannot guarantee the convexity of our objective function, the algorithm is only guaranteed to converge to a local optimum [13] . As a result, we advocate initializing the algorithm using multiple restarts from different initial starting values [14] . We find in our empirical studies that numerically optimizing the log-likelihood using the Nelder-Mead algorithm is stable and can repeatedly deliver a consistent optimal solution from several distinct initializations.

With a robust method for computing the MLE, we similarly can obtain interval estimates numerically. Recall for large samples, asymptotic normality of the maximum likelihood estimator can be used toward constructing approximate confidence intervals [15] , [16] . Denotingθ ML θ ML θ ML as the maximum likelihood estimate for θ θ θ, its asymptotic distribution is given byθ 

We now assess the proposed methodology via simulation. The experiment aims to recover the true parameters used to generate the synthetic data from the immigration-death model given data at only a set of discrete observation times.

We design the simulations to resemble the real data in one of our case studies, with the ground truth parameters and observation schedule being chosen so that the shape and scale of the simulated outbreaks reflects that of the OC hospitalization data we seek to analyze.

The daily number of hospital and ICU beds occupied in the county-level OC data is obtained by aggregating the individual numbers of hospital and ICU beds occupied as reported by each major hospital in the county. We replicate this in our simulations by first generating the number of hospital and ICU beds occupied for several hospitals, and then summing across hospitals to obtain population-level hospitalization data that is comparable to the OC data we will study. As it is possible that some hospitals do not report hospitalization data on a particular day, we additionally extend our simulation study to assess the robustness and potential bias of the proposed framework in the presence of underreporting. Note that we might expect the simulated hospitalization data for each hospital in the absence of any underreporting to resemble the UCI Health hospitalization data, while the aggregated hospitalization data in the presence of underreporting should resemble the noisier hospitalization data at the county level.

Simulation Procedure To establish notation, letH tj andĨ tj be the number of hospital beds and ICU beds occupied on the t th day in the j th hospital and [0, 1, 2, . . . T max ] be the discrete times at which we observe the process with [0, T max ] being the total observation period. We first generate the daily number of hospital beds occupied for each of the m hospitals. To reflect scale of the population-level data in our case study, we simulate individual hospital occupancies so that the total number of occupied beds across hospitals matches our dataset. To this end, we begin by samplingH 0j from a multinomial distribution such that m j=1H 0j = H O (0), the total number of hospital beds at the start of our observation period in the OC data. For each day in the observation period, we next sample the difference

The daily number of hospital beds occupied,H(t), during the entire observation period for m hospitals is then obtained by cumulatively summing d tj .

Conditional on the simulated hospital bed countsH tj , we generate the daily numbers of ICU beds occupied from an immigration-death process with rates

We use the Gillespie algorithm [17] to generate realizations of the process defined by the rates in Equation (8) After independently simulating 1000 instances ofH tj , we simulateĨ tj conditional on eachH tj as described above.

We then aggregate each of the 1000 instances ofH tj andĨ tj ,

to obtain the daily total hospital and ICU beds occupied respectively.

In order to capture possible underreporting by hospitals, we next consider a misspecified simulation by randomly sampling a number γ t ∼ U nif {0, 1, 2} of hospitals to be "dropped" on each day in the observation period. These dropped hospitals are chosen uniformly at random from the 25 hospitals, and are omitted from the sums (9) before aggregating the countsH tj andĨ tj .

Results For both studies with or without underreporting, we calculate the transition probabilities using Equation (7) with N = 2 11 , and numerically optimize the log-likelihood from ten random restarts under the Nelder-Mead algorithm. Relative Error Figure 1 : Violin plots of relative errors in estimated immigration and death rates from the simulation study. The corresponding coverage rates of the 95% confidence intervals are mentioned above the violin plots.

range, results and uncertainty estimates especially pertaining to µ must be interpreted conservatively when only coarse population level data are available.

We now analyze hospitalization data provided by UCI Health, the University of California Irvine's hospital system, collected from August 02, 2020 to December 31, 2020. UCI Health provides us access to anonymized patient-level data which includes the duration of ICU stay for each patient that was diagnosed with COVID-19 and admitted at UCI Health. The availability of such "line-list" is a rare exception, and in this case study will serve as a ground truth to validate our method, allowing us to check against inference from aggregating the individualized data. We explore a number of candidate models with variables of interest including the number of hospital beds and ICU beds occupied by patients who have a confirmed diagnosis for COVID-19 and the daily cumulative test positivity rate. Model 3, with

was found to be the best model fit for the UCI Health data via AIC. Figure 3 appearing in Appendix C.

Since the COVID-19 pandemic was rapidly evolving in OC and the plot of daily total number of hospital and ICU beds occupied by patients showed distinct trends over the entire observation period, we decided to divide the data into After dividing the data into three phases, we apply our immigration-death framework to estimate all parameters, with a focus on interpreting the per-patient ICU stay as the reciprocal of the death rate µ. We first consider analysis under a time homogeneous model, where we assume that model parameters mentioned in Table 1 are constant across all three phases i.e., from March 29, 2020 to November 15, 2020 and fit all models based on this assumption. Next, we fit the parametrized models described in Table 1 to each phase separately, referring to this case as the time inhomogeneous model.

The results from fitting each model to the data from OC in the time homogeneous setting suggests that model 3, with

is preferable in terms of achieving the has the lowest AIC on the county-level data. 

attains the lowest AIC value during phase 1, while model 3 with

attains the lowest AIC values during phases 2 and 3. Further details of the AIC values under candidate models are reported in Table 3 appearing in Appendix D.

During phase 1, our primary parameter of interest, average per-patient ICU stay, was estimated to be 2.56 days, with 95% confidence interval [1.92, 3.45]. Our inference suggests this average stay increased to 5.99 [4.17, 10 .0] days during phase 2, and then slightly decreased back to 5.00 [3.33, 10 .0] days during phase 3. This is a novel finding from fitting the population level time series, and while there is no validation data such as line-list records available for the county-level study, our validation on the UCI Health data gives us assurance in our results. Though they are not all as directly interpretable, estimates and confidence intervals for the remaining parameters across all three phases are depicted in Figure 2 . We see that phase 2 reflected the largest uncertainty in estimates, while the majority of estimates are statistically significant. A detailed tabulated summary of these results is reported in Table 2 appearing in Appendix D.

In order to further justify model fit, we perform a validation study using model-based simulation from the estimated parameters of the time inhomogeneous model. We also assess the division of the total time period into three phases by repeating the validation study using the estimated parameters for the time-homogeneous model. For each setting, µ ,β H ,β P respectively, for the best selected model during the three phases in OC.

we ran 100 repeat simulations of the number of ICU beds occupied conditional on the hospital bed counts. We save the 2.5 th and 97.5 th percentile of the simulated daily number of ICU beds occupied to create pointwise middle 95% quantiles of the predictive distribution.

Figures 3a, 3b, and 3c display these 95% quantiles for the predictive distribution of the simulated daily ICU beds for the best time homogeneous model for OC, the time inhomogeneous models for OC, and the homogeneous model for UCI Health, respectively. We see that in the time-inhomogeneous case, the predictive interval contains the observed number of beds occupied through the majority of the time series plot. This suggests that the model successfully balances flexibility with model complexity. On the other hand, this is in contrast to the time homogeneous model which leads to a much more pronounced discrepancy between the predictive intervals and the observed counts.

We propose a novel likelihood-based framework to conduct inference on discretely observed data from an immigrationdeath process that can be adapted to a flexible class of covariate-dependent rates. By deriving the likelihood function for partially observed count data, we are able to numerically perform maximum likelihood estimation and obtain confidence intervals. Moreover, because this framework gives access to criteria such as AIC and BIC, it further allows us to perform model selection among a class of parametric models. By applying our framework, we are able to estimate time-varying parameters describing the per-patient immigration rates into ICUs as well as interpretable durations of ICU stay from population level data collected in OC. Moreover, a similar application restricted to data from UCI Health offers validation compared to direct estimates from available line-list data.

In our case study, Phase 1 coincides with the rise in initial infections in California, a time when there was no standard hospital policy for tackling the rise in infections. As the virus had never been encountered before, hospitals may have erred on the side of caution, generous about sending patients to the ICU regardless of the extent of their conditions. As a result, patients who were being sent to the ICU may have been recovering and and getting discharged at a faster rate, thus reducing the average per-patient intensive care occupancy.

By the beginning of the second phase, hospitals were better adapted to handle COVID-19 patients. As a result of enhanced guidelines from the state government, hospitals may have become more selective in admitting patients into the ICU. Perhaps those in more vulnerable groups, such as elderly and patients with underlying medical conditions that increased their probability of developing severe symptoms and complications, were the main groups being admitted into the ICU, leading to an increase in the average stay. The average stay in the ICU decreased slightly during phase 3, which directly follows the summer peak. In general, due to some level of underreporting in the data, it is safer to treat our results as slight underestimates of the average per-patient ICU stay. Nonetheless, our validation study using model-based simulation is encouraging, and future work may explore observation models under various emission distributions to directly model the effect of under-reporting.

In models developed to forecast hospital occupancy during the COVID-19 pandemic, it is common to obtain parameter estimates using knowledge from other regions. Unlike other mechanistic models which may calibrate model parameters based on outside information, this paper proposes methodology for estimating the average per-patient ICU stay using only hospitalization data. Together with the mechanistic model's ability to better capture the underlying ICU stay dynamics by allowing the rates to be functions of covariates, this suggests that the proposed inferential framework can be embedded as extended compartments within models such SEIR [18] . Hospitalization and ICU compartments have been key additions within disease models, both because hospitalization data is an observable feature of a pandemic and because of the importance of forecasting hospital system occupancy during an outbreak ( [19] , [20] ). An advantage of our immigration-death framework is that it makes use of surveillance data instead of requiring patient-level data, which is typically difficult to obtain and requires HIPPA compliance. In our current case studies, our model selects daily positivity rate and total number of hospital beds occupied as relevant covariates. Future work may investigate the impact of additional covariates such as the introduction of vaccines and therapeutics, or the implementation of more efficient hospital protocols on ICU admittance rates. It would also be fruitful to consider non-parametric estimation of changes in the death rate, and to examine its effect on the accuracy of the estimated per-patient ICU stay.

[20] Olga Morozova, Zehang Richard Li, and Forrest W. Crawford. One year of modeling and forecasting COVID-19 transmission to support policymakers in connecticut. Sci Rep, 11:20271, 2021.

The data and open-source code to reproduce all results and experiments described in the article are made available at the authors' webpages, while the data for the case study of Orange County is publicly available 5 .

A Deriving the closed form solution for φ i (s, t).

We can differentiate the PGF to calculate the Kendall's partial differential equation for an immigration-death process ([9] ). On differentiating both sides of equation (4) in the main text we get,

Thus, the Kendall's partial differential equation for this process is,

The solution of equation 10 gives an explicit expression for the PGF. The solution of equation 10 is,

B Algorithm to simulate the number of ICU beds occupied for a set of hospitals over a fixed observation period. Figure C1 : Immigration-death model of ICU capacity. If k beds are occupied in the ICU, then after a random waiting time there can be either k +1 or k −1 beds occupied. The immigration-death model is characterized by an immigration rate α(t) that does not depend on the current beds occupied in the ICU and a death rate µ k = µk that is proportional to the current beds occupied in the ICU. 

A model to forecast regional demand for COVID-19 related hospital beds. medRxiv

Inference for nonlinear dynamical systems

Transition probabilities for general birth-death processes with applications in ecology, genetics, and evolution

Inference for reaction networks using the linear noise approximation

Birth/birth-death processes and their computable transition probabilities with biological applications

Parameter estimation for multivariate population processes: a saddlepoint approach

A linear noise approximation for stochastic epidemic models fit to partially observed incidence counts

A First Course in Stochastic Processes

Applied Probability

Calculation of the equilibrium distribution for a deleterious gene by the finite Fourier transform

Likelihood-based inference for discretely observed birth-death-shift processes, with applications to evolution of mobile genetic elements

A simplex method for function minimization

Numerical Optimization

Algorithms for Optimization

The consistency and ultimate distribution of optimum statistics

Exact stochastic simulation of coupled chemical reactions

Likelihood-based inference for partially observed stochastic epidemics with individual heterogeneity

Forecasting hospitalization and ICU rates of the COVID-19 outbreak: An efficient SEIR model. Bull World Health Organ, E-pub

This work was partially supported by NSF grants DMS 2030355 and DMS 1936833, as well as a Donald Bren School of Information and Computer Sciences Exploration Award.