Hospital Bed Occupancy Sensor

Forecasting local hospital bed demand for COVID-19 using

On-Request Simulations

Angelo D’Ambrosio1, Raisa Kociurzynski1, Alexis Papathanassopoulos1, Fabian

B¨Urkin1, Hajo Grundmann1, And Tjibbe Donker∗1

Germany.

Abstract

For hospitals, realistic forecasting of bed demand during impending epidemics of infectious diseases is essential to avoid being overwhelmed by a potential sudden increase in the number of admitted patients. Short-term forecasting can aid hospitals in adjusting their planning and freeing up beds in time.

We created an easy-to-use online on-request tool based on local data to forecast COVID-19 bed demand for individual hospitals. The tool is flexible and adaptable to different settings, and it is based on a stochastic compartmental model for esti- mating the epidemic dynamics and coupled with an exponential smoothing model for forecasting.

The models are written in R and Julia and implemented as an R-shiny dash- board. The model is parameterized using COVID-19 incidence, vaccination, and bed occupancy data at customizable geographical resolutions, loaded from official online sources or uploaded manually. Users can select their hospital’s catchment area and adjust the number of COVID-19 occupied beds at the start of the simulation.

The tool provides short-term forecasts of disease incidence and past and forecasted estimation of the epidemic reproductive number at the chosen geographical level. These quantities are then used to estimate the bed occupancy in both general wards and intensive care unit beds. The platform has proven efficient, providing results within seconds while coping with many concurrent users.

By providing ad-hoc, local data informed forecasts, this platform allows decision- makers to evaluate realistic scenarios for allocating scarce resources, such as ICU beds, at various geographic levels.

Keywords— Infectious diseases, COVID-19, Dashboard, On-request forecast-

Introduction

The speed of community transmission of virulent pathogens during epidemics or pandemics of infectious diseases has the potential to overwhelm healthcare systems . Hospital beds are usually in limited supply, because most hospitals run at near- full capacity under normal circumstances. A sudden increase in patients needing hospital care might thus not be met with the required available beds if the growth in demand is unexpected. Reliable forecasting of bed demand during epidemics is therefore essential for hospitals to free up beds in time.

However, epidemic trajectories can differ considerably regionally , and because patients are usually admitted to hospitals in their own region , the bed demand can widely vary between hospitals as well. Consequently, nationwide forecasts of bed demand may not reflect the individual hospital’s situation, and may fail to predict overload problems. We, therefore, developed a bed demand forecasting model based on local epidemic trajectories in the hospital’s individual catchment area and its current bed occupancy.

To be widely applicable, this model needs to be able to produce bespoke fore- casts for each hospital in a country for which basic levels of data are available. This capability requires a flexible forecasting system because of the vast number of po- tential combinations of parameter choices for all individual hospitals, in particular concerning the possible choices for catchment areas per hospital, the number of cur- rently occupied beds, and patients’ lengths of stay on various wards. Such a variety of parameters enables the extreme forecast flexibility required to produce ad-hoc reports and forecasts for all situations of specific hospitals.

We made the use of this model accessible for non-technical users by creating an online interactive dashboard. This platform allows individual users (i.e., hospital managers) to enter the data and parameters related to their specific context and requirements and produce ad-hoc bed occupancy estimates and predictions. End user-specific forecasting of bed demand is a unique concept, as on-request forecasting is often avoided because of the computational challenges of multiple concurrent user modelling platforms.

Here we show how some of these computational challenges can be solved. Careful streamlining of the model and its technical implementation can help produce epi- demic forecasts in a reasonable time, allowing users to explore the potential future epidemic trajectory without being restricted to the developer’s viewpoint.

We will first discuss the forecast model itself, its data and parameter require- ments, as well as the different modules it contains. Each module is designed to guide the user through the steps and assumptions needed to produce the forecast.

Then we will discuss the technical implementation of the online dashboard and the technologies that ensure usability for multiple concurrent users.

Finally, We Will

present performance statistics of the dashboard collected between December 2021 and March 2022. The code used to implement the dashboard, in the version used to prepare this manuscript, is available at: https://github.com/QUPI-IUK/Bed-demand-forecast/releases/tag/v.0.5.6.

General Model Structure

The complete forecast model consists of 5 main modules: 1) Data loading and nowcasting, 2) Reproduction number estimates, 3) Vaccination coverage forecasting, 4) Incidence model, and 5) Care path model. The calculations done in the first two modules rely only on the loaded data and user input, while modules 3-5 also build further on the results of each previous module (see figure 1).

This Hierarchical

structure implies that any change to parameters in a module updates the outputs in all dependent modules. However, by guiding the user through the modules in sequential order, the number of required computations can be drastically reduced.

User Input

Figure 1: Model structure showing each of the modules, the data carried over between modules, and the possible user input in each module.

Data Structure And Definitions

The model needs data on infection incidence, vaccine doses, bed occupancy in general wards and intensive care units, and a simulation window for which a forecast is required. Although we developed this model for hospitals in Germany, based on data provided by the Robert Koch Institute and the Deutsche Interdisziplin¨are Vereinigung f¨ur Intensiv- und Notfallmedizin (DIVI) , the model can be applied to any country that has the following data available (preferably if in remote-accessible,

2. Administered Vaccinations, Vd,G(T)

3. bed occupancy on general wards BGW,g(t) and ICU BICU,g(t) where t denotes the time in days since a reference date before the start of the pandemic, g denotes a geographical subdivision (in the case of Germany, this is a district, i.e., Landkreis or Stadtkreis), and d denotes the dose of vaccine administered (1st, 2nd, or 3rd/booster). Incidence and vaccination data are needed for each day

3

of the period or reference, while merely the value at the start date of the forecast is needed for bed occupancy. The user defines the catchment area of interest, selecting geographical areas from the available list (Kall), into the catchment set (Ku). After this, the input data is

G∈Ku

BICU,g(t). The start of the simulation is denoted as TS, for which the incidence, number of vaccinations, and bed occupancy are taken from the (last data-point of the) data, the first forecasted day is thus TS + 1. The length of the forecast is defined as L; the last forecasted date is therefore given as TE = TS + L. The size of any class of individuals in the model is given in absolute numbers of individuals; thus S(0) = N, with N being the population size.

Because many of the parameters in the dashboard’s models are continuously dis- tributed (as Exponential, Gamma, or Weibull), and the model uses discrete time steps, we need to discretise the parameter distributions. For each continuous distri- bution, with probability density function (PDF) fc(t), the discretised PDF is defined

T

fc(x)dx. User input for all distributions is given in terms of the mean µ and standard deviation σ, and distribution parameters are then determined by moment matching.

For a gamma distribution, this means that f(x, α, β) = f(x, α = µ2/σ2, β = µ/σ2), for a Weibull distribution, f(x, k, λ) = f(x, k = ( σ

Μ)−1.086, Λ = Μ/Γ(1 + 1/K)), And

for an exponential distribution f(x, λ) = f(x, 1/µ).

Vaccination Model

The vaccination model consists of two parts.

In The First, We Forecast The Num-

ber of people vaccinated over the length of the forecast for each of the doses (1st, 2nd, and 3rd/booster dose), while in the second part we convert these administered vaccinations into population-level protection against transmission.

The first doses (V1(t)) are assumed to be administered at the same rate as during

V1(I)/7

The second dose is assumed to be administered at a fixed time delay after the first dose. This time delay ∆TV is extracted from the data as the maximum time for

4

holds. The second doses are then a direct reflection of the number of first doses

∆Tv Days Ago,

V2(t) = V1(t −∆TV ). The daily administered third/booster doses are assumed to be a continuation of the mean daily third doses over the previous week:

I=Ts−6

V3(i)/7. This is in line with the forecasting of the first doses, with the exception that the cumulative number of third doses is not allowed to exceed the cumulative number of second doses a certain delay (∆TB) ago, such that

V2(I)

has to hold.

Population Protection

Each dose is assumed to have a time-dependent additive effect on the population- wide protection against infection. At the individual level, Gd(τ) denotes the propor- tion of individuals no longer susceptible to infection τ days after the administration of dose d.

This function thus serves to simulate the delay between vaccine ad- ministration and maximum vaccine protection at the individual level. We assume no waning of immunity, such that Gd(∞) = Ed, with Ed the vaccine effectiveness against transmission elicited by doses d.

The additive effect on vaccine effectiveness of each additional dose is defined as

Ed

∗= Ed−Ed−1; 1. consequently, the additive increase in an individual’s protection

Ed

. Gd(τ) is input as the cumulative distribution function (CDF) of a normal dis- tribution with mean µ(G1) = 15, µ(G2) = 15, µ(G3) = 7, and standard deviation σ(G1) = 3.8, σ(G2) = 6.5, σ(G3) = 3.8 as default values.

The population protection at time t is then given by,



. which is used in both the R forecasting step and the incidence model.

R Estimation And Forecasting

In order to inform the transmission process part of the incidence model, we need to forecast the reproduction number of the pathogen over the model timeframe. This forecast is based on the observed time-varying effective reproduction number (Re(t)), and implemented either as the static continuation of the most recent observation, or forecasted using an ETS (Error, Trend, Seasonal) model based on the last 100 days of observations, depending on the user’s preference.

Furthermore, We

provide the option to incorporate the effect of a Variant of Concern (VoC) on the development of the R0(t). In all cases, we use estimated values of Re(t) from the available data up to the simulation start date TS, Re(t ≤TS), and forecasted values after the simulation start date (Re(t > TS)).

We estimate the time-varying effective reproduction number Re(t) using the Cori et al. R estimation method as implemented in the EpiEstim R package (version 2.2-4), with a given non-parametric Serial Interval (SI). This SI distribution denotes the number of days between equal disease events (e.g., onset of symptoms) of directly connected cases. The default SI is assumed to be Gamma distributed with a mean of 5 days and a standard deviation of 4.9. These values were chosen to reflect estimates of mean SI from multiple studies with a relatively high standard deviation.

By keeping track of the cumulative number of infected individuals, as well as the previously calculated population protection (GP (t)), we can calculate R0(t) from Re(t), given that Re is related to R0 through the susceptible population:

S(T)

.

Ets Model

To produce a dynamic forecast of R0(t), we used an ETS model which underlies an exponential smoothing method . In general, exponential smoothing methods generate point forecasts of time series using weighted averages.

For an observed time series y1, y2, ..., yn the point forecast for h periods ahead is yt+h. When the forecast is based on all data up to time t, it is denoted as y∗ t+h|t.

For a method with an additive trend, the point forecast results from the level of the

T+H|T = Lt + Hbt

which gives the forecast equation. While lt and bt are based on the level and trend (slope) at t −1 and are given

Lt = Αyt + (1 −Α)(Lt−1 + Bt−1),

bt = β∗(lt −lt−1) + (1 −β∗)bt−1.

6

This approach is also known as Holt’s linear method.

Values For The Initial

states l0 and b0 as well as the smoothing parameters α and β∗are estimated from the observed data, with 0 < α, β∗< 1. Increasing α gives more weight to the more recent observations and less weight to observations that lie longer ago. In our case, α and β∗were set to 0.25 and 0.15, respectively.

While the exponential smoothing methods generate point forecasts, the underly- ing statistical models, in addition, generate prediction intervals. To generate such an ETS forecast model, we need to specify the probability distribution for the residual at time t, et. For additive errors, we assume the residuals are normally and indepen- dently distributed with mean 0 and variance σ2, notated as et = εt ∼NID(0, σ2).

Expressing the error as εt = yt −lt−1 −bt−1 ∼NID(0, σ2) we can substitute εt into

Lt = Lt−1 + Bt−1 + Αεt,

bt = bt−1 + βεt. with β = αβ∗. These expressions are referred to as the state-space models as they include an equation for the observation and equations for the unobserved states (level, trend). Each observation yt in the ETS model therefore includes a component for the level (or smoothed value), an error component (E), a trend component (T), and seasonality (S). In our case, we use the ETS(A,A,N) model, which has additive

Errors, Additive Trend And No Seasonality;

We use a log(R0(t)) time series of a 100 days, i.e. log(R0(TS −99...TS)), as input for the ETS model. The logarithm serves to forecast the relative changes in R0 and avoids negative R0 forecasts. Per model iteration, we pick a random quantile for the prediction interval and use this as the trajectory for R0(t > TS).

Variant Of Concern Implementation

To model the impact of new pathogen variants emerging in the population , we need to consider how a variant of concern (VoC) has an advantage over the back- ground viral population. We considered two possible mechanisms: 1) through higher transmissibility, increasing the base R0(t), and 2) by immune evasion, effectively in- creasing the size of the susceptible class (S(t)).

In case of increased transmissibility, we increase the basic reproduction number of the background viral population R0,B(t) by the relative added fitness advantage

A,

R0,V (t) = (1 + A)R0,B(t). Note that the previously calculated R0(t) is the average of the background and VoC reproduction numbers weighted by the proportion of cases caused by either

R0(T) = Ρb(T)R0,B(T) + Ρv (T)R0,V (T),

= (1 −ρV (t))R0,B(t) + ρV (t)(1 + A)R0,B(t). We can calculate R0,B(t), used as the general reference reproduction number in

7

or we can directly calculate this from the observed Re(t) value

Re(T)N

((1 −ρV (t)) + ρV (t)(1 + A)) S(t). The proportion of cases caused by the VoC (ρV (t)) over time can be described with a sigmoidal function. A similar approach to estimating fitness advantages for new variants has been used by others .

This Is Done By Assuming Both The

background population and the VoC have a static growth rate unaltered by added immunity through new cases (i.e., the relation between their growth rates is fixed over time). We can calculate these growth rates (r) as a function of R0, using the

Σ

. This means that the number of infections from both background variants and the

Σ

t.

T

. By defining the total number of infected individuals as IT (0) = IV (0) + IB(0), and

Σ

t . This shows that we only need the relative added fitness advantage (A) in combi- nation with a starting proportion of the variant (ρV (t = tρ)) at a certain reference date (tρ) to describe the proportion VoC over time. In the dashboard code, the reference date for the proportion of VoC does not overlap with the starting date of the incidence, vaccination, and bed occupancy data. In other words, with I(tI = 0) and ρV (tρ = 0), tI̸ = tρ, because ρV (tI = 0) = 0, making a precise estimate of ρV (t)

8

impossible, or at least impractical. We, therefore, set tρ to a later time point at which the ρV is identifiable. The added relative fitness advantage of the VoC can then be estimated from data about the proportion VoC among all isolates per week in the given catchment area, available through for instance genomic surveillance efforts. The user can then set the reference date, starting proportion, and added relative fitness advantage based on such analyses to be used in the forecast.

Immune Evasion

In case the VoC has an advantage over the background variant through (partial) immune evasion, the number of individuals susceptible to the VoC (SV (t)) would be greater than for the background variant (SB(t)), giving it a higher Re(t) even if the R0(t) would remain the same.

We define immune evasion (ϵ) as the proportion of individuals immune against infection with the background variant that are not immune against infection with the VoC. Thus, the number of individuals immune to infection with the VoC scales linearly with the number of individuals immune to infection with the background variant, if very few or no infections with the VoC have occured yet.

To correctly assess the added fitness advantage at the VoC reference point (tρ), we need to calculate the increase in the number of individuals susceptible to infection with the VoC relative to those susceptible to the background variant. For simplicity, we assume that the cumulative number of infections with the VoC is negligible relative to the cumulative number of infections with the background variants at tρ, that vaccination is aimed at the background variants, and immunity against the VoC is only determined by cross-immunity (1−ϵ) from immunity against the background variants.

Sv (Tρ = 0) = N −(1 −Ε)(N −Sb(Tρ = 0))

Given that Re,B(t) = R0(t)SB(t)/N and Re,V (t) = R0(t)(1+A)SV (t)/N, the rel- ative added advantage by immune evasion, given the pre-existing immunity against

= (1 + A)Sv (Tρ)

SB(tρ) −1. As A∗is also the observed fitness advantage, the user will need to make an assumption about what part of this advantage is caused by immune evasion or by increased transmissibility on the reference date. Although the reference date (tρ) and the start of the simulation (TS) are not necessarily the same, since the user can define tρ independently of TS, from here on, we only consider the case where TS = tρ, for clarity and brevity.

Incidence Model

The incidence forecast model is based on a stochastic implementation of a Susceptible- Infected-Recovered (SIR) model. The model describes the number of new cases as a function of the population still susceptible (S(t)) to infection and of infection pres- sure (P(t)) exerted by all those currently infectious, and the total population size (N). The size of the components is given as absolute numbers of individuals; thus S(0) = N. The infection process is governed by the forecasted basic reproduction number R0(t ≥TS) computed in the preceding modules of the framework.

If no VoC is defined, the model describes the dynamics of a single variant (the background variant). Conversely, the model converts into a two-strain model if a VoC is defined, keeping track of the infected individuals for either variant as well as the total number of individuals susceptible to either variant.

Single Variant Model

The size of the susceptible class is determined by both the cumulative number of infected individuals and the population level protection through vaccination, as de-

!

. The infection pressure exerted on S(t) is then given as the weighted number of infected individuals in the preceding serial interval,

H(I)I(T −I)

These are thus all past cases that are able to cause new cases on the current day t, weighted by the probability that they do so i days after they were reported themselves. Given the current R0(t) and P(t), the expected mean total number of newly infected individuals on day t+1 is I(t + 1) = R0(t)P(t).

The transmission rate per infected individual (The probability of infecting each of the other individuals in the population in a unit of time), is given as β(t) = R0(t)/δN. Since we use infection pressure P(t) calculated over the entire SI, we can set δ to 1, such that β(t) = R0(t)/N.

The infection pressure exerted on each individual within the population is then

Pi(T) = 1 −(1 −Β(T))P (T),

which for low numbers of infected individuals (and therefore low P(t)) can be sim- plified to PI(t) = βP(t). The newly infected individuals (I(t+1)) are then randomly picked from a binomial distribution I(t + 1) = Binom(p = PI(t), N = S(t)).

Voc Implementation

If a VoC is defined, we need to track the size of both the susceptible and the infected class of individuals for both the background variant and the VoC. This is first done retrospectively, for t ≤TS, to determine the size of the classes at TS.

10

With an estimate of the VoC proportion over time (ρV (t ≤TS)), the past infec-

Iv (T ≤Ts) = ⌊Ρv (T)I(T)⌉,

where ⌊...⌉denotes rounding to the nearest integer. The number of infections with the background variants thus follows as IB(t) = I(t) −IV (t). This also gives us the exact size of the susceptible class for both the background variant and the VoC

!

. Note that in the absence of new cases, the susceptible classes still need to be calculated forward in time before the incidence model is actually run, as forward vaccinations reduce their size.

The forward, prospective, model of the infection process, for t > TS, is then performed for both the background variants and VoC in the same way as with a single variant, splitting I(t), P(t), PI(t), and S(t) in a background variant (∗B) and VoC (∗V ) class. The final epicurve then combines both variants again, I(t) = IB(t) + IV (t).

Care Path (Within-Hospital) Model

To forecast bed occupancy in the hospital, the main endpoint of this work, we first simulate the path of each single patient admitted to the hospital through the general ward, ICU, and step-down units (care path model), as described in Donker et al (). Then we join this model with the predicted incidence to forecast the hospital admission rate and the bed occupancy.

We use the following parameters to model the care path of individual patients through the hospital. • Length of stay (LoS) distribution on general wards (GW)

• Los Distribution On Step-Down Units (Sdu)

• proportion of patients being transferred from GW to ICU • proportion of patients being transferred from ICU to SDU It is possible to add the proportion of patients directly admitted to ICU (HI) to this list, but this is usually estimated from the data (see below).

Once all admissions and discharges are determined, we track the total number of patients on each ward type on each day of the forecast, presenting the final result of the forecast. Note that, despite modelling the general ward and step-down unit separately, we combine both wards into a single number of beds occupied on the general ward in any other step, because the occupied beds in both cases are usually reported as part of the same type (i.e., occupied non-ICU beds).

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.