Design of COVID-19 staged alert systems to ensure healthcare capacity with minimal closures

Community mitigation strategies to combat COVID-19, ranging from healthy hygiene to shelter-in-place orders, exact substantial socioeconomic costs. Judicious implementation and relaxation of restrictions amplify their public health benefits while reducing costs. We derive optimal strategies for toggling between mitigation stages using daily COVID-19 hospital admissions. With public compliance, the policy triggers ensure adequate intensive care unit capacity with high probability while minimizing the duration of strict mitigation measures. In comparison, we show that other sensible COVID-19 staging policies, including France’s ICU-based thresholds and a widely adopted indicator for reopening schools and businesses, require overly restrictive measures or trigger strict stages too late to avert catastrophic surges. As proof-of-concept, we describe the optimization and maintenance of the staged alert system that has guided COVID-19 policy in a large US city (Austin, Texas) since May 2020. As cities worldwide face future pandemic waves, our findings provide a robust strategy for tracking COVID-19 hospital admissions as an early indicator of hospital surges and enacting staged measures to ensure integrity of the health system, safety of the health workforce, and public confidence.


Supplementary
: Diagram of the epidemiological model with compartments, transitions, and rates, with definitions provided in the enumerated notation that follows. The model has ten copies of such diagrams for each age-risk group pair, which interact to determine the (unmarked) rate of transition from compartment S to E. Such transitions are governed by the transmission rate β but also contact matrices and the current risk stage.
We assume that the threshold for each stage is constant over time and is denoted by i for stage i. The goal is to 48 determine level i so that the total expected cost of our trigger policy is minimized, while respecting epidemiological 49 dynamics and hospital capacity using a probabilistic constraint. We approximate the expectation and probabilistic 50 constraint by using Monte Carlo simulation to form a set of scenarios, indexed by ω, and by using binary variables to 51 count the number of paths that violate capacity and to determine the stage for each sample path at each point in time. a ∈ A set of age groups {0-4y, 5-17y, 18-49y, 50-64y, 65y+} r ∈ R risk groups {low, high} i ∈ I predefined stages {1 (red), 2 (orange), 3 (yellow), 4 (blue)} ω ∈ Ω set of scenarios

Parameters
Epidemiological parameters: β unmitigated transmission rate σ rate at which exposed individuals become infectious τ proportion of exposed individuals who become symptomatic ρ A rate at which pre-asymptomatic individuals become asymptomatic ρ Y rate at which pre-symptomatic individuals become symptomatic γ A recovery rate from asymptomatic compartment γ Y recovery rate from symptomatic compartment γ a H recovery rate from hospitalized compartment for age group a γ a ICU recovery rate from ICU compartment for age group a P proportion of pre-symptomatic transmission Y HR a,r percent of symptomatic infectious that go to the hospital for age-risk group a, r η H hospitalization rate after symptom onset ω A infectiousness of individuals in IA relative to IY ω a,r P P 1−P τ (Y HR a,r /η H +(1−Y HR a,r )/γ Y )+(1−τ )ω A /γ A τ /ρ Y +(1−τ )ω A /ρ A : infectiousness of pre-symptomatic individuals relative to IY for age-risk group a, r π a,r γ Y ·Y HR a,r [η H −(η H −γ Y )Y HR a,r ] : rate-adjusted proportion of symptomatic individuals who go to the hospital for age-risk group a, r p IH percent of patients directly going to the general ward of the hospital HICU R percent of general ward patients who get transferred to ICU η a ICU ICU admission rate after admission to the general ward for age group a ν a H γ a H ·HICU R [η a ICU −(η a ICU −γ a H )HICU R] : rate-adjusted proportion of general ward patients transferred to ICU for age group a µ a rate from ICU to death for age group a ICU F R a percent of hospitalized that die for age group a ν a ICU γ a ICU ·ICU F R a [µ a −(µ a −γ a ICU )ICU F R a ] : ICU fatality rate-adjusted proportion for age group a φ a ,r ,a,r i,t expected number of daily contacts from (a , r ) to (a, r) at time t under stage i N a,r population of age-risk group a, r i the stage-specific threshold (level) for the daily hospitalization rate X i,t,ω 1 if the system is in stage i at time t for scenario ω; 0 otherwise Y i,t,ω 1 if the trigger statistic lies in range of stage i at time t for scenario ω; 0 otherwise V i,t,ω 1 if Y i,t,ω = 1 prior to spending 14 days in the previous stage for scenario ω; 0 otherwise Z ω 1 if healthcare capacity is exceeded in scenario ω; 0 otherwise We refer to Supplementary Table 10 for further details on model parameters. We first define the epidemiological 54 transition dynamics in the following equations for all ω ∈ Ω. These dynamics largely follow the formulation used  For simplicity, we write the finite-difference Eqs. [1] in a deterministic form. They become stochastic, and require indexing by ω, because binomial random variables replace terms like σE a,r t,ω ; here the binomial random variable has parameter n = E a,r t,ω and σ serves as the "success" probability. This construct is pervasive throughout right-hand side terms in Eqs. [1]. In addition to these "micro" stochastics there are "macro" stochastics because we model σ, ω A , γ A , and γ Y as random variables that are subject to a Monte Carlo draw at time 0 of the simulation. We discuss this further in Supplementary Method 3. replacing ICU a,r t,ω and B ICU by ICU a,r t,ω + IH a,r t,ω and B.

90
We have formulated the model with daily time periods, which simplifies notation for computing the seven-day 91 moving average of new admissions, the logic behind the rule of spending two weeks in a stage, etc. That said, in 92 implementation we use ten time steps per day, which suffices for the fidelity of the epidemiological dynamics in 93 Eqs. [1]. 94 We generate 300 scenarios for Ω according to the procedure described in Supplementary Method 3. After obtaining 95 a set of optimal triggers, we generate another 300 scenarios, Ω , to assess performance. Those two sets of scenarios 96 are similar to "training" and "testing" data used in statistics and machine learning. As a result, it is possible to obtain 97 thresholds that meet the probabilistic constraint under Ω, but violate that constraint when tested using Ω , although in 98 our experience such violations are both rare and modest.

99
Model [3] is a large-scale stochastic mixed-integer nonlinear program. Problems of the scale we consider cannot currently be solved using commercial integer programming software. We approximately solve the model using a grid search procedure. For a fixed set of thresholds, i , i ∈ I, which satisfy inequalities [3m], we can run the simulation model for all ω ∈ Ω in parallel, applying transition dynamics [ ∀ω ∈ Ω. [3o]

100
The cost coefficients C 1 , C 2 , C 3 , and C 4 are important parameters in the optimization model [3]. The foremost driver 101 in the probabilistic constraint requires that we respect ICU capacity. Policies that satisfy that constraint specify a set of 102 allowable thresholds. We choose C 1 C 2 C 3 C 4 so that the resulting policy is "hierarchical." Effectively, the 103 model first tries to find policies that minimize expected time in stage 1 and satisfy the probabilistic constraint on ICU 104 capacity. Among those it tries to minimize the expected time spent in stage 2, etc. For our specific problem instances, 105 we first computed results with C 1 = 10 3 , C 2 = 10 2 , C 3 = 10 1 , and C 4 = 10 0 . Then we solved the model again 106 with C 1 = 10 6 , C 2 = 10 4 , C 3 = 10 2 , and C 4 = 10 0 . The obtained optimal thresholds and the detailed results are 107 identical in both cases and correspond to what we report in Table 2 in the main text. This suggests that the factor of 10 108 for C i /C i+1 is sufficiently large to induce the desired hierarchical property. This overall cost structure represents the 109 preference of avoiding increasingly strident alert levels, and has been proven effective in the City of Austin in practice, 110 as discussed in the main manuscript and as further detailed in Supplementary Discussion 1.  We define four baseline contact matrices, H, S, W, and O, to describe the contact frequency between age groups otherwise. [4] The indicator 1 {off day} takes value 1 if the day is a weekend or holiday and is otherwise 0, and a similar indicator 128 accounts for school closures. When a high-risk group, along with those 65 years and older, is involved either on the 129 "giving" or "receiving" end of a contact, Eq.
[4] assumes reduced transmission via the cocooning coefficient, c i .

130
The following are key dates during the pandemic in Texas, and some define time blocks, which we use in estimating Travis County (which includes Austin) banned gatherings of more than 100 people [11,12].

143
• July 17, 2020: Time point in hospitalization data suggesting a change in dynamics.

144
• August 9, 2020: The last day of observed data from the hospital system used in estimating changes in ICU 145 dynamics.

146
• August 20, 2020: First day students returned to residence halls at the University of Texas at Austin.

147
• October 7, 2020: The last day of observed data used in estimating model parameters. 148 We assume that there are six time blocks denoted by T j for j ∈ {1, 2, 3, 4, 5, 6} as defined in Supplementary Table 3.

149
They guide fitting of transmission-reduction parameters, κ and c, and certain dynamics in use of the ICU and hospital 150 duration, as detailed below.  We model the hospitalization dynamics, including proportions of hospitalized requiring the ICU and durations in the general ward and ICU, using data from a multi-facility hospital system serving the central Texas region, including Austin, Texas ("Hospital system data"). While we model differences based on five age groups, we assume the same hospital dynamics in different hospital systems after a patient is admitted across Austin and Houston due to similar medical standards. Conditional on being admitted to the hospital, we observe a decreasing trend in the probability a patient is admitted to the ICU throughout the time horizon, which holds for both direct admissions to the ICU and patients who are first admitted to the general ward. Among patients who enter the general ward and are then admitted to the ICU, their duration of stay in the general ward, determined by η ICU , grows over time. For each time block, T j , we assume a constant η ICU,j and further assume a constant daily decrease, r j , on both of the fractions, p IH and HICU R: along with a similar decrement across boundaries of the blocks. We use duration times for each time block from the 152 hospital system data to estimate η a ICU,j and fit r j , with the estimated parameters in Supplementary Table 4.  Using the hospital system data, and consistent with the transition diagram in Supplementary Figure 1, we define the ICU duration for a patient as the time between their admission to the ICU and their discharge from the hospital. The reality is more complex as ICU patients typically return to the general ward prior to discharge from the hospital, and iterations between the two units, driven by a patient's health status, can also occur. Therefore, the reported duration in the ICU leads to over estimating ICU utilization and under-estimating that of the general ward. To handle this in our model, we introduce two constant parameters, α ICU and α H , to better estimate durations in the ICU and general ward and better represent their respective utilization: where γ 0 H , γ 0 ICU , and µ 0 are obtained from the hospital system data, with each row corresponding to an age group in ascending order: with units of day −1 .

154
The bulk of the epidemiological and hospitalization parameters are specified above or are detailed in Supplementary Tables 10 and 11, with the latter obtained from the literature or information collected from local healthcare 156 agencies. The time blocks are specified in Supplementary Table 3. Given these, we estimate 14 parameters, but with 7 157 degrees of freedom, as we detail below. We perform the fit of the deterministic SEIR model in Eqs. [ Houston MSA. By minimizing a weighted sum of least-square errors, we estimateκ j andĉ j , j = 1, 2, . . . , 6, α H , and 161 α ICU , using SciPy/Python [13] via scipy.optimize.least_squares. We use the "hat" notation on κ and c 162 to distinguish, for example,κ j for time block j (see Supplementary Table 3) from κ i , which corresponds to stage i per 163 Eq.
[4], and show the mapping shortly. 164 We minimize where IH t , ICU t , and H t denote the estimated IH t , ICU t , and H t obtained through Eqs. [1]; w ICU and w H are 165 scaling constants; and the sum is over t ∈ T 1 ∪ · · · ∪ T 6 . We assume w ICU = 1.50 and w H = 7.58, as those 166 values approximate magnitudes relative to that of the general ward. To obtain a parsimonious model, we useĉ 1 = 0, 167ĉ 2 =ĉ 3 =κ 2 ,ĉ 4 =ĉ 5 =κ 4 =κ 5 andĉ 6 =κ 6 , which reduces the number of estimated parameters from 14 to 7. The 168 rationale is that there was effectively no cocooning during the initial time block T 1 , and thus we setĉ 1 = 0. Because  The physical distancing parameters for each stage are mapped to κ i and c i for i ∈ I based on the historical implementation of the policy. We set: In Eq.
[6a] we set κ 1 (red) to the strictest observed level of transmission reduction, June 26-August 19, 2020.  evitably, some stochastic sample paths will yields hospitalizations that diverge from observed data, including sample 193 paths in which spread quickly terminates after initializing with a single infectious person.

194
Here we describe how we select scenarios, sampled from the marginal distributions of σ, ω A , γ A , and γ Y , and from the micro-stochastics process, in order to generate a set of 300 scenarios, indexed by ω ∈ Ω and used in Eqs. [1]. Because the projections at earlier time points guide the trigger optimization model at later time points, the goal is to start the latter with scenarios, ω ∈ Ω, that are consistent with observed hospitalizations up to that point in time. In order to evaluate the quality of scenario ω in this sense, we use where IH is the mean of the IH t values. Algorithm 1 summarizes the sampling procedure to generate |Ω| = 300 195 scenarios.

196
Algorithm 1 uses the six time blocks T j , j = 1, 2, 3, 4, 5, 6, given in Supplementary Table 3. We start with 197 simulating a sample path ω during the first time block T 1 (see line 4-5), and compute the corresponding R 2 value at 198 the end of T 1 . If the scenario, ω, has the desired quality over T 1 , i.e., R 2 ≥ 0.85, we continue simulating the path ω.  For the purpose of visualizing results we use spaghetti plots that detail 300 sample paths, and we also use a single 224 "representative" path. A simple least-squares selection of a path from the collection of 300 paths arguably does not 225 yield a representative path. Because the timing of peaks in paths differs, the simple least-squares path is "flatter" 226 through the cloud than most paths. Here we describe how we select a representative path from the collection of 300. 227 We define a selection criterion based on metrics of interest, such as the total number of people hospitalized in the general ward over the duration, the peak number of people in the hospital, and the time of that peak. We let IH t,ω , ICU t,ω , and H t,ω denote the estimates of IH t , ICU t , and H t for each sample path, ω ∈ Ω, and we first obtain the squared deviations from the observed data, defining Z ω,obs as follows: where w IH and w ICU are scaling constants. Here, t ranges over all days up to October 7, 2020, for which we have 228 observed data. We then define metrics to represent the statistical properties of the scenarios as in Supplementary 229 Table 6. The first and second row represent the set of scenarios for Austin and Houston, respectively. Rather than showing the rates, we show 1/σ, 1/γA, and 1/γY , which correspond to durations in days. We sample from the marginal distributions of each of the four parameters, and Algorithm 1 accepts sample paths consistent with hospitalizations, while also sampling a binomial number of transitions between compartments in the SEIR-style model. While there are modest deviations from the nominal triangular distributions for durations (1/σ, 1/γA, and 1/γY ), the distribution of the relative infectiousness of an asymptomatic individual (ωA) deviates significantly from its nominal distribution. Here, σ is the rate at which exposed individuals become infectious; ωA is the infectiousness of asymptomatic relative to symptomatic individuals; and γA and γY are recovery rates from the asymptomatic and symptomatic compartments.
Analogs of the standardized metrics for those hospitalized in the general ward (IH) from Supplementary Table 6 231 are also computed for those hospitalized in the ICU (ICU ), and new daily hospital admissions (H) to form a total of 232 12 such metrics. We select the representative scenario ω such that 233 ω ∈ arg min ω∈Ω w std 12 Z ω,std + w obs Z ω,obs where w std and w obs are positive weights that sum to one. This method is used to select the representative paths 234 Supplementary Figure 3: Scatter plot of the random variables 1/σ, ωA, 1/γA, and 1/γY for the scenario set Ω selected by Algorithm 1. The left and right scatter plots represent the set of scenarios, Ω, for Austin and Houston, respectively. We sample independently from the marginal distributions of each of the four parameters, and the algorithm accepts sample paths consistent with hospitalizations. This can induce dependencies as shown, for example, between 1/γY and ωA. Here, σ is the rate at which exposed individuals become infectious; ωA is the infectiousness of asymptomatic relative to symptomatic individuals; and γA and γY are recovery rates from the asymptomatic and symptomatic compartments.  Two key factors play a role the Monte Carlo acceptance-rejection sampling procedure of Algorithm 1. One is the 238 nature of the stochastics, i.e., macro-stochastic and micro-stochastics. The second is the initial "burn-in" period of 239 February 28-March 23, 2020 (for Austin), which plays an important role, in part because we start with a single 240 infectious individual on February 28th. In our analysis, we generate two sets of scenarios that are similar to "training" 241 and "testing" data used in statistics and machine learning. We use 300 scenarios for each set, and for Austin this 242 required generating 918,220 total scenarios; i.e., the acceptance rate is about 0.07%. Of those scenarios, all but In Section 2 of the main text we compare four alternative policies to our optimized four-stage trigger policy. First, 256 we test a benchmark which is again optimized and differs from the four-stage policy only in that it has access to just 257 two stages of mitigation, red and yellow. In this case, we require a single threshold to identify when to toggle between 258 the two stages. In our second benchmark, we again optimize in the same way as the four-stage policy, but we do so 259 using a 0.95-level probabilistic constraint to respect total hospital capacity rather than ICU capacity. 260 We assess two additional benchmarks. Instead of triggering based on the seven-day moving average of daily hos- The two left-most columns indicate the relative weight on transmission reduction factors for the orange and yellow stages. So, the last row corresponds to yellow, while the 50-50% row corresponds to equal weight on the orange and yellow values. The two right-most columns show the resulting transmission reduction coefficients for the low-risk population (κ) and the high-risk population (c) for Austin.

271
We assume the ACS can be set up within two days; i.e., the ACS will be ready to use two days after the trigger days of unmet demand, under these two scenarios, respectively. We would expect demand for ACS beds to last up to 297 two months, and to decrease with higher levels of adherence.

298
Sensitivity Analysis

299
We report additional results on the sensitivity associated with public compliance with recommended physical 300 distancing, or lack thereof. Supplementary Figures 8a and 8b   From left to right, the three optimal policies are derived to prevent overwhelming ICU demand using either a four-level or two-level alert system or to prevent overwhelming inpatient demand using a four-level system. As benchmarks, we evaluate policies implemented in France [14] and proposed as gating criteria for relaxing measures and opening schools in [15]. For Houston, the orange and red thresholds for the Percent ICU policy translate to 300 and 600 COVID-19 ICU cases, respectively; the yellow, orange, and red thresholds for the Incidence policy translate to 710, 7100, and 17750 new cases, respectively, assuming that one in ten cases is reported. We implemented each policy in our stochastic SEIR model fit to hospitalization data for the Houston, Texas MSA assuming the reported COVID-19 ICU and inpatient capacities of 1000 and 4500 beds, respectively. Outcomes are based on 300 stochastic simulations of COVID-19 transmission and healthcare burden from October 7, 2020 through September 30, 2021 under each policy.     Table 9: Optimized ACS policies under four adherence scenarios for the Austin area. The scenarios assume that the public do not fully comply with enacted alert levels, resulting in a constant transmission rate lying somewhere between those expected under the orange and yellow stages. The percents in the second column indicate the relative weighting of the orange and yellow rates; for example, the moderate-high rate is obtained via a 75-25 weighting of the two rates, respectively. The ACS policies track the moving seven-day average of COVID-19 hospital admissions and assume that the ACS will be ready to accept patients two days after triggering. The launch trigger is the derived threshold for opening the ACS; the target capacity is the estimated total ACS beds needed, based on the criteria that 95% of projections remain under ACS capacity; probability launched is the fraction of simulations in which admissions exceed the launch trigger; the 95th percentiles of days used and peak demand provide upper bounds on the projected demand. Coupled optimization of the the trigger and capacity in discrete steps of 100 beds and 10 daily admissions, respectively, along with the resulting slack between peak demand and target capacity can produce small non-monotonicities in optimal triggers, as seen in the middle two rows. We assume the ACS remains open for 100 days once triggered, which is more than sufficient to cover the estimated 95th percentile of days required (penultimate column    The reduced ICU capacity (from 331 to 200 beds) led to a significantly lower (stricter) threshold for moving from the 397 orange to red stage (a seven-day average COVID-19 hospital admission of 50 versus 120 ω P : infectiousness of individuals in pre-symptomatic and preasymptomatic compartments, relative to symptomatic and asymptomatic compartments  austin-travis-county-alternate-care-site-ready-take-patients.