A hospital demand and capacity intervention approach for COVID-19 in the UK

Following CVAal20 , in this section we present a simple data-driven susceptible-exposed-infected-recovered-dead (SEIR-D) model as demonstrated in Figure 1. The model breaks down the typical infectious compartment of an SEIR-D model into two compartments, one strand to model individuals who become infectious with COVID-19 and will be going to hospital ($I$), and the other strand to model individuals who will not need to go to hospital and thus remain undetected by hospital healthcare requirements ($U$). Splitting the infectious compartment this way allows us to describe the number of individuals who are in hospital with COVID-19 at any time, which allows us to fit the hospital data to the model to be able to infer information about the number of infectious cases in the community. This approach allows us to circumvent the issue of parameter identifiability and estimation due to the lack of obtainable and usable data regarding those who are infectious and asymptomatic.

In CVAal20 we used specific regional datasets to calibrate the model. As part of the national COVID response, all the National Health Service (NHS) hospitals in England treating COVID-19 patients submitted a Daily Situation Report (SITREP) to NHS England. The data associated to Sussex NHS Trusts were extracted and combined to weekly death counts from the Office for National Statistics (ONS) with COVID-19 reported as the underlying cause of death, which we then received and fitted. To be precise, by Sussex we mean the collective term for geographies pertaining to the counties of East Sussex and West Sussex in South East England. The SITREP contained counts of daily admissions, daily discharges, and the beds occupied daily, whilst the death dataset contained the number of deaths recorded outside of hospitals and the number of deaths recorded within hospitals.

Whilst the ONS death data are publicly available, the SITREP data in general are not, but, given the national need for data, the UK government produced the Coronavirus Dashboard. It provides users with the ability to look at different metrics of COVID-19, such as hospitalisations and deaths, for different regions and provides an API for users to download the datasets themselves. The granularity of the data depends on the size of the region the data are required for, typically all data are available for each nation of the UK, but hospitalisations and deaths are split depending on their geographical location. Indeed, hospitalisations are recorded using NHS regions and NHS trusts, and deaths are recorded using local authority regions. The unfortunate difference between using the Coronavirus dashboard and the SITREP is that the Coronavirus dashboard does not contain the number of daily discharges, and it does not differentiate place of death like the ONS weekly death registry does. In order to apply the approach presented in CVAal20 we adapted the Coronavirus dashboard data in the following way: we calculated the proportion of hospital COVID-19 deaths to the total number of COVID-19 deaths for the region being considered from the ONS weekly death registry and applied this proportion to the deaths dataset acquired from the Coronavirus Dashboard to give us a proxy dataset on deaths in hospital and deaths outside of hospital. Next, we then used the deaths in hospital with the admissions and beds occupied to find a proxy dataset for the discharges, since the offset of beds occupied between each day depends on the number of admissions, discharges and deaths that day. We note that, although the proxy death datasets might closely track the ONS weekly death registry, the discharges dataset is very noisy due to noise accumulating from multiple sources. For this reason, we decided to apply a 7-day rolling average to the discharges. We note here that the Coronavirus dashboard has access to patients in mechanical ventilation beds and so we could include a compartment that describes the high-dependency unit (HDU). We decided not to do this since we do not have access to the number of patients who have died in HDU, we only have access to the number of patients who died in hospital, and so parameter identifiability would be an issue.

For illustrative purposes, we only consider the geographical regions of North West England, South East England and the nation of England. We chose these as both NHS region and local authority region are the same.

The mathematical model is a data-driven generalisation of a simple temporal epidemiological dynamic system of ordinary differential equations, governed by the interactions depicted in Figure 1, supported by non-negative initial conditions

Here, the dot above the notation denotes the time derivative. In this setting, $N$ denotes the total regional population. $S(t)$ denotes the proportion of the total population $N$ who are susceptible to the disease, COVID-19. Susceptible individuals become exposed to the disease, i.e. they are carrying the disease but are not currently infectious, to form the $E(t)$ subpopulation at rate $\lambda(t)$ which represents the average infectivity. The rate $\lambda(t)$ is the product between $\beta$, the average transmission rate, and the probability of meeting an infectious person $(U(t)+I(t))N^{-1}$. The $E(t)$ subpopulation is an incubation compartment and further evolves in two ways. A proportion of $E(t)$ becomes infectious but, in the spirit of the hospital healthcare demand, remains undetected with probability $p$, denoted $U(t)$, at a rate $\gamma_{E}$, or a proportion of $E(t)$ becomes infectious and will require hospitalisation in the future with a probability of $1-p$, denoted $I(t)$, at a rate $\gamma_{E}$. The subpopulation that does not require hospitalisation can either progress to recover with a probability of $1-m_{U}$, at rate $\gamma_{U}$, to form the recovered population, denoted by $R_{U}(t)$, or die with a probability of $m_{U}$, at rate $\gamma_{U}$, to form the dead population, denoted by $D_{U}(t)$. Considering the infectious population that will be going to hospital, these individuals will become hospitalised, denoted by $H(t)$, and thus be in hospital care at rate $\gamma_{I}$. We assume that once a patient has been admitted into hospital, they no longer transmit to other non-COVID-19 patients, visitors or workers. This could be due to appropriate isolation and use of personal hygiene measures in hospital, however we do note that transmission in hospital is not unheard of. For the model, this means that we assume transmission is only conducted in the community. Once in hospital, patients can evolve in two separate pathways, a proportion of the hospitalised population can fully recover at rate $\gamma_{H}$ to form the subpopulation denoted by $R_{H}(t)$. Alternatively, if they can not recover, then they die while in hospital at rate $\mu_{H}$, to form the dead population denoted by $D_{H}(t)$. We note that $\gamma_{H}$ and $\mu_{H}$ are closely related to the length of stay in hospitals.

As is standard for epidemiological models of this nature, $\beta$ denotes the average transmission rate, $\gamma_{E}^{-1}$ denotes the average incubation time, $p$ denotes the average proportion of infectious individuals who will not require hospital treatment, $\gamma_{U}^{-1}$ denotes the average infectious period for those not needing hospital treatment, $m_{U}$ denotes the infected fatality ratio for undetected cases, $\gamma_{I}^{-1}$ denotes the average infectious period from becoming infectious to being admitted to hospital, $\gamma_{H}^{-1}$ denotes the average hospitalisation period for those who recover and $\mu_{H}^{-1}$ represents the average hospitalisation period for those who die. For this model, using the method of next generation matrices DHR10 , we derive the formula for the basic reproduction number $\mathcal{R}_{0}$ as the following

Using this, we calculate the effective reproduction number $\mathcal{R}_{t}$ as

The effective reproduction number $\mathcal{R}_{t}$ is often referred to as the $R$ number. We briefly explain the fitting procedure presented in CVAal20 , we refer the interested reader to CVAal20 for further details. First, we utilise the linear relationship between the model description of hospital discharges and hospital deaths and use linear regression analysis to calculate the ratio of discharges to deaths. One can see that, in terms of the model and its parameters, the daily discharges can be written as

and the daily hospital deaths can be written as

The linear regression allows us to estimate the ratio between $\gamma_{H}$ and $\mu_{H}$. Next, we rewrite equations (2.1)–(2.5) in terms of the data, we call this the “observational” model, whereby (2.6)–(2.9) are not considered since they are cumulative representations of the compartments in (2.1)–(2.5). The observational model is formed of compartments that are available from the data. Indeed, one can see that, in terms of the model and its parameters, the daily admissions can be written as

the daily deaths outside of hospital is

and daily discharges as in (2.12) above. Thus, the observational model, as presented in CVAal20 , is

We solve the observational model, compute (2.12), (2.14) and (2.15) and compare against the datasets given. This allows us to use a maximum likelihood estimation approach by minimising a negative log-likelihood with some constraints on the initial conditions and on the effective reproduction number. These constraints are needed to keep the feasible region of the parameters close to the realistic sets. Without these constraints, the peak of the datasets could be explained by being close to herd immunity (i.e., a lot of infections have already occurred before the lockdown) or a large value of $\mathcal{R}_{t}$ combined with a small value of $p$. We note here that (2.13) does not need to be used in the observational model due to the linear regression and the fact that the resulting log-likelihood functions from the regression and observational model are independent. This means that there is one less parameter to infer from solving the observational model.

Since we are considering the data from the first day of lockdown, we also need to infer the initial conditions. We have not conducted a formal investigation into the resulting log-likelihood, but it is clear that there is a continuous dependence between initial conditions and the parameters, see CWVDM21 for a comprehensive discussion. In practice, we see this manifest as an issue to fit $p$ - if the first guess of initial conditions and parameters is not close to the “true” values, then it is $p$ which changes in value the most. Although in reality, $p$ is characterised by how COVID-19 affects individuals, from demographics such as age, ethnicity, gender, however we speculate this value should not change drastically between regions. In view of this, in the parameter estimation process we fix $p$ to be the same throughout the estimation for the three regions.

From the model and using the fitting procedure described above, we gain the parameters in Table 1 and initial conditions in Table 2, with a demonstration of the fit for beds occupied demonstrated in Figure 3. The equivalent figure for the Sussex region fit can be found in CVAal20 in Figure 2. We note that the infected fatality ratio and average hospitalisation period for those who recover are similar across all regions, but the average hospitalisation period for those who die and the value of $\mathcal{R}_{t}$ when lockdown commenced are quite different. The varied values of $\mathcal{R}_{t}$ could be explained by the amount of infections in each region, it was reported at the time that Sussex and the South East escaped quite lightly on the number of infections, whilst the North of England did not. This is shown in Figure 2, whereby we have taken the model output of the beds occupied and divided it by the population size for each region and multiplied it by 100. Proportionally, the North West saw almost double the amount of patients as the South East. The higher value of $\mathcal{R}_{t}$ for North West can also be seen by looking at the gradient of the beds occupied. There could be a myriad of reasons for this, such as the demographic of the population or the geography of the region. One should also note that, as we came out of the first lockdown and went into the first iteration of the tiered system, it was cities in the north which first started to show signs of resurgence (such as Manchester). Using $\gamma_{H}$ and $\mu_{H}$ for each region, one can calculate that the estimate of the probability of discharge for England, South East and North West was 58.03%, 61.57%, 55.31% respectively. As for the reason why, one may speculate that this is due to the number of patients in the hospitals being higher in the north which puts pressure on the health system.

As a reference to compare how the interventions are working, we use this section to demonstrate the “do-nothing” approach, which is simply to let the disease take its course. This will provide us with statistics to compare to the interventions later on to demonstrate their effectiveness. In reality, we are aware that this approach will not be implored in practice and an intervention will occur, as has happened all over the world with national level lockdowns and social distancing measures. Since the parameters presented in Table 1 depict a lockdown scenario, we scale $\beta$ appropriately to several values to establish an epidemic, i.e. so that $\mathcal{R}_{0}>1$. An increase in average transmission rate can simply be interpreted as more individuals meeting each other and spreading the disease. In Table 3 we measure the maximum hospital capacity needed as a percentage of the total population, the day in the simulation that maximum is reached and the percentage of dead individuals at the end of the outbreak, for each region. As intuitively expected, as $\mathcal{R}_{0}$ increases the percentage of maximum number of patients in the hospital increases, the day of that peak is sooner and the percentage of dead individuals increases. Interestingly, the day of the peak does not seem to change even though the hospitalisation parameters are all quite varied, we suspect that this is somehow related to the initial conditions and the fact that $\mathcal{R}_{0}$ is not effected by the parameters that describe hospitalisations. We see that the percentage of maximum beds occupied is worse in the South East, followed by England and then the North West respectively. This is to be expected since the average hospitalisation periods (for deaths) are longer in the respective order which implies more people in hospital at any one moment. However, the percentage of dead individuals is higher in the North West, followed by England and then the South East respectively. This is also to be expected since the probability of death in hospitals is larger in the respective order whilst the mortality probability outside of hospitals are very similar for each of the regions. In Figure 4 we demonstrate the effective reproduction number $\mathcal{R}_{t}$ of the simulation using the England parameters. This shows us that, when $\mathcal{R}_{0}$ is larger, the actual outbreak is much shorter in length and reaches much smaller values of $\mathcal{R}_{t}$. This description follows the notion that the larger the value of $\mathcal{R}_{0}$, the more aggressive the disease is following the exponential growth of those who are infectious, as can be seen in Figure 4 by the steep decline in $\mathcal{R}_{t}$. In Figure 5 we demonstrate a comparison of the percentage of beds occupied for each of the regions and for some values of $\mathcal{R}_{0}$. We note that in Figure 5 we truncate the simulation to make the visualisation of the hospitalisations easier.

In the following sections we want to investigate how one can use hospital capacity as a measurement for whether interventions are put into place. We aim to model the situation where an intervention is triggered when hospital capacity is almost full, and then finish the intervention when the hospital capacity has reached an “opening” threshold of significantly lower patients. For ease of demonstration we will simply set a threshold that once breached will trigger the intervention, in some cases this will mean that the capacity is “breached” - we deal with this in subsection 3.5. We denote the state of being in an intervention using the notation $\ell$, namely if $\ell=1$ then we are in an intervention, otherwise we set $\ell=0$. Using this, we describe the average transmission rate as

where $\beta$ is the average transmission rate associated to the first lockdown deduced from the parameter inference for each region, $C_{\mathcal{R}_{0}}$ is a scaling constant to give the initial value of $\mathcal{R}_{0}$ wanted, and $[\cdot]$ are the Iverson brackets K92 defined as

We describe the intervention in the following recursive way

which is to say that the intervention is triggered at time $t$ when the number of patients in hospital goes above an upper limit $H_{u}$, and the intervention stays triggered until the number of patients in hospital goes below a lower limit $H_{l}$. Initially $\ell$ is set equal to 0. The values of $H_{u}$ and $H_{l}$ are regionally dependent since they depend on the maximum capacity of all the hospitals in a region. Given the maximum hospital capacities in Table 3, we fix $H_{u}:=0.0012N$, to guarantee at least one lockdown for each $\mathcal{R}_{0}$ chosen, and vary the opening up threshold between $H_{l}:=0.025H_{u}$ and $H_{l}:=0.5H_{u}$. We will present the results from all the regions in various tables, however graphically we will restrict the figures to the South East region so the figures are not a visual burden. To this end, we present Figure 6 which is the only comparison of the interventions for each region. The simulations all look reasonably similar and follow a similar trend, the noticeable differences come from when one region goes into an intervention (light grey coloured line) whilst the others don’t. This explains some of the interesting numbers in Tables 4–5. We go into more detail in the following sections.

There are some caveats to this study we want to highlight. It is somewhat unrealistic that transmission would immediately revert to normal amounts after an intervention MB20 , however we do not consider this here. Another aspect to consider with this study is that we are considering the capacity of all the hospitals in a region, due to modelling constraints and data access. This means that we are assuming hospitals can move patients throughout each region due the physical constraints of each hospital, in response to the bed capacity of each individual hospital.

In this section, we will vary $\mathcal{R}_{0}$ whilst fixing $H_{l}:=0.25H_{u}$ or $H_{l}:=0.5H_{u}$. In Table 4 we measure the percentage of dead individuals for each region for varying values of $\mathcal{R}_{0}$. In Figure 7 we demonstrate the percentage of beds occupied in the South East region whilst in Figure 8 we demonstrate the effective reproduction number over the simulation. The light grey colour in Figure 7 represents the time when the regions are in an intervention, the same time period is represented by the gaps in the lines in Figure 8.

Comparing Table 3 to Table 4 we can see that the intervention has dramatically decreased the percentage of total deaths for larger values of $\mathcal{R}_{0}$, as expected. By looking at Figure 7 we can see that as we increase $\mathcal{R}_{0}$, the number of interventions needed increases and also the length of the initial intervention increases. This is due to the number of future patients in the exposed compartment $E$ and the infectious compartment $I$. This is emphasised by the fact that the initial intervention is sooner for a larger value of $\mathcal{R}_{0}$ due to the increased average transmission. It can also be seen that the time between each intervention decreases as $\mathcal{R}_{0}$ increases. We can also see, by comparing Figure 4 with Figure 8, that the epidemic actually lasts significantly longer with interventions included. These observations are realistically expected, however there are some results which are not necessarily expected or intuitive such as the percentage of dead individuals in the South East region decreasing from $\mathcal{R}_{0}=1.5$ to $\mathcal{R}_{0}=1.6$ when considering $H_{l}=0.25H_{u}$ but not decreasing for when $H_{l}=0.5H_{u}$. This is not specific to this value of $\mathcal{R}_{0}$, rather to the circumstance that this simulation finds itself after the final intervention. Namely, at this stage of the simulation for $\mathcal{R}_{0}=1.6$, herd immunity has almost been reached, i.e. $\mathcal{R}_{t}$ is only slightly larger than 1. A value of $\mathcal{R}_{t}$ slightly larger than 1 means that although there is still an increase in the number of infectious individuals, the rate of that increase is much slower comparatively to an $\mathcal{R}_{t}$ value of, say, 1.5. At this stage, the number of incubating and infectious individuals in the case of $\mathcal{R}_{0}=1.6$ are small which means that the final bump in the simulation for $\mathcal{R}_{t}=1.5$ is much larger than the respective bump for $\mathcal{R}_{t}=1.6$, as can be seen in Figure 7. One can see a similar behaviour happening between $\mathcal{R}_{0}=1.4$ and $\mathcal{R}_{0}=1.5$ when considering $H_{l}=0.5H_{u}$. This interplay between parameter values, herd immunity and $\mathcal{R}_{t}$ is difficult to analyse and demonstrates that intuition is not necessarily enough when forecasting.

In this section, we will be experimenting by changing $H_{l}$, the lower threshold of patients that signals the ending of the intervention, and fixing $\mathcal{R}_{0}=1.5$ or $\mathcal{R}_{0}=1.6$. In Table 5 we measure the number of interventions needed, the length of each intervention (measured in days), the day of the initiation of each intervention and the percentage of the total deaths at the end of the epidemic. We demonstrate a few of the simulations in Figure 9 and we demonstrate the effective reproduction number over each outbreak in Figure 10.

Comparing Table 3 to Table 5 we can see that again the intervention has dramatically decreased the percentage of total deaths compared to the “do-nothing” approach. Interestingly, changing the threshold $H_{l}$ mostly does not have that much of an effect on the percentage of total deaths, unlike the difference we saw in Table 4 when changing $\mathcal{R}_{0}$, but it does have a large effect on the length of the outbreak. The effect of $\mathcal{R}_{t}$ is slightly different in this section, namely the sudden drop in the percentage of deaths is not caused primarily by a well-timed re-opening. Instead, the jump in percentage is due to the fact that an extra intervention has been used. Considering the South East region, and focusing on $\mathcal{R}_{0}=1.5$, we can see that the percentage of deaths drops from 1.1682% using $H_{l}=0.1875H_{u}$ to 1.000% using $H_{l}=0.5H_{u}$, and by Figure 9 or 10 that there is one intervention using $H_{l}=0.1875H_{u}$ whilst there is two interventions using $H_{l}=0.5H_{u}$.

The outcomes of our study so far imply that the timing and the lengths of interventions are extremely important. Getting closer to herd immunity when ending an intervention has the potential to save a huge number of lives. However, calculating $\mathcal{R}_{t}$ in real life is in general very challenging which leaves the process of timing for herd immunity difficult. One also notices that, although the percentage of total deaths decreases with an intervention, the length of most of the interventions is large due to the criteria set. This is mainly due to the fact that the average hospitalisation period is large and that the scenario we are simulating means that intervention will be in place until hospitals go from full capacity to between 2.5% and 75% capacity. Fortunately, we see that as the target capacity percentage increases, the percentage of total deaths does not increase dramatically, and the length of interventions decreases from the best part of 4 months to under 1 month. Similarly, as the outbreak progresses, one would expect the average hospitalisation length to decrease, since awareness of the disease and treatment gets better, as well as an increase in resources and the development of vaccines. This final point is important as it means realistic interventions can be implemented as circuit breakers and still maintain a large decrease in the number of total deaths. However, one aspect of this which is overlooked in this study is the potential for nosocomial outbreaks, whereby the probability of an outbreak increases with a larger number of infectious patients.

In this section we calculate what the highest capacity threshold $H_{u}$ is, for different values of $\mathcal{R}_{0}$, so that the hospitals do not go over their maximum capacity $H_{max}$. In particular, this approach can be used as an early warning system by outlining when interventions need to be enforced to maintain a manageable capacity. This approach can complement scenario-based forecasting approaches by giving a range of indicators of when to open up further capacity in hospitals or introduce an intervention which can be tracked against with incoming data daily. Ultimately, in practice, healthcare systems will want to utilise their capacity appropriately whilst not impacting the general public with an intervention, and so optimising the difference between the resource capacity and the maximum number of patients per parameter set is important.

In this section, we only consider the situation of finding $H_{u}$ using one intervention, that is to say (3.3) becomes

with (3.2) the same. We only need to consider one intervention here because in most situations, once the hospitals have opened back up again, there will be another spike in capacity, however one can reason that this second spike is essentially the same as the first spike, just with slightly different parameters. For a specified value of $\mathcal{R}_{0}$, we look to find the root of the following function

where we note that $H$ depends on $H_{u}$. Since $H$ is continuous, this function is continuous and will always have a root provided $\mathcal{R}_{0}$ is chosen high enough such that the maximum of $H$ without an intervention is at least as large as $H_{max}$. Figure 11 depicts the situation where $\mathcal{R}_{0}$ is chosen appropriately for the South East region. For this section, we take $H_{max}:=0.0012N$.

We demonstrate the results in Figure 12. This figure should be read in the following way: if we know $\mathcal{R}_{0}$, then we look to determine that an intervention should be initiated when the healthcare system reaches a certain percentage of the maximum resource capacity. We can see that as $\mathcal{R}_{0}$ increases, the cap on the occupancy for the intervention decreases, which is to be expected. This measurement is useful, as, depending on the value of $\mathcal{R}_{0}$, one can calculate what the cap is on occupancy and thus can judge when to take steps to initiate an intervention.

In this section we will present the agent-based version of the equations in (2.1)–(2.9). In particular, we look to draw parallels between the two schematics in Figures 1 and 13. In a similar manner to Figure 1, Figure 13 describes the way an agent moves between states and should be read as a flowchart, where we have labelled the states similarly to demonstrate some of the parallels. In Figure 13 the circles represent the state an individual can be, and the diamonds represent a decision. The arrows from a states to a decisions are an event in the epidemic, and the arrows from a decision to a state or decision are the outcome of the decision. Some arrows have an outcome of a decision with a parameter in parentheses next to them, these parameters are probabilities of traversing a different section of the flowchart and represent the same parameter as when there is a fork in states in Figure 1, such as the probability of not going to hospital $p$. Some arrows have an event with a state in brackets next to them, these are to distinguish between the different rate parameters as the states represent similar characteristics, such as the multiple types of being infectious $U$ and $I$.

For ease of exposition, we consider the time unit of the simulation to be days. Each event decision is decided by drawing one pseudo-random uniform number from the interval $[0,1]$, let us call this value $x$, and comparing it to the probability of the event on a given day, let us call this value $d$. We can describe this process as the following categorical function

We note that for each event, a new $x$ is generated, and each $d$ is different. The events can be split into two categories, transmission and progression. We begin by describing the transmission event, event: contact. For ease of exposition, each day we assume that all agents make a fixed $C$ random contacts with a probability of a successful transmission $a$, where $\beta=Ca$. For a susceptible agent, of those $C$ contacts, we denote $\kappa$ to be the number of contacts they make with an infectious agent, namely an agent in state $U$ or $I$. Then, we have that

where $X\sim Geo(a)$ and $X$ has support on $\{1,2,\dots\}$. The remaining events are all progression events, whereby the length of time an agent has spent in a state is important. We will describe “event: incubating”, but note that the same description can be used for each of the events using their respective rate parameters. Let $l$ represent the number of days an agent has spent in the $E$ state, then we have

where $Y_{E}\sim Exp(\gamma_{E})$. One can see that event: infectious $[U]$ gives $Y_{U}\sim Exp(\gamma_{U})$, event: infectious $[I]$ gives $Y_{I}\sim Exp(\gamma_{I})$, and event: hospitalisation gives $Y_{H}\sim Exp(\gamma_{H}+\mu_{H})$. We note that this is not the standard methodology for agent-based models, typically the length of stay for each agent is not stochastically generated, but instead each day an agent runs a Bernoulli trial using the probability associated to the event of progression. In our framework, it is obvious to see that compartments are exponentially distributed, and there is an obvious way to change that assumption if needed without having to describe the underlying mathematics differently. For example, if it was known that the incubation event was in fact akin to an Erlang distribution rather than an Exponential, then we can update the agent-based approach quite simply by updating the random variable $Y_{E}$, however in the equation-based approach we would need to add more compartments (the additional number depending on the distribution’s shape parameter) CDE18 .

For the remainder of this study, we will focus only on the outcomes using the South East parameters. In order to compare the agent-based approach to the standard deterministic SEIR-D approach, we run 100 Monte Carlo simulations and take the average number of agents in each state for each day. We note that typically 100 Monte Carlo simulations are not enough but we see very little variation in the average from using as low as 5 iterations. Due to the stochastic and numerical discretisation, we adjust the rate parameters by adding a correction term in the form of

where $\Delta t$ denotes the time interval being considered in ratio of days, e.g. $\Delta t=0.5$ would be considering half days, see Appendix A for justification. We note that $\Delta t$ needs to be chosen appropriately so that the denominator in $c$ is non-zero, however this is a reasonable assumption since typically $\gamma<1$ (as it is a rate and its reciprocal in normally larger than one) and $\Delta t\leq 1$ since we typically consider simulating each day or several simulations per day. Using the parameters defined in Table 1, with the associated fitted initial conditions, taking $\Delta t=0.25$ results in Figure 14, whereby we see a clear agreement between the two approaches and the data. We note that this result is also due to the fact that $N$ is large, and the initial conditions are suitably large, which combined are often called the mean-field assumption (or the thermodynamic limit). From here onwards, we will fix $\Delta t=1$ for speed but we note in Appendix A that using $\Delta t=1$ and the correction term results in outputs sufficiently close to the SEIR-D model.

In this section, we look to apply the techniques and parameter changes in Section 3 to the agent-based approach. We know that $\beta$ is multiplicatively formed of the average probability of transmission $a$, and the average number of contacts per day $C$. When changing the transmission rate, $\beta$, we will consider two different approaches. We look to investigate whether just changing the probability of transmission $a$ or average number of contacts per day $C$ have different impacts on the outcome of the epidemic. One can reason that a change in probability could be due to new strains of an infectious disease or a change in public behaviour, such as wearing masks or increased personal hygiene. Changing the number of contacts could be due to a national level intervention, such as lockdowns or closing schools. We note that we whilst it would be scientifically ideal to keep $a$ fixed, due to the integer nature of contacts $a$ has to be adjusted to maintain the same $\mathcal{R}_{t}$ value across the measurements, thus making sure the tables are comparable.

By using the “do-nothing” approach, we see that in fact changing either parameter does not make an overall difference to the outputs, which is to be expected given the mean-field limit of the agent-based approach to the compartmental model. From here onwards, we will fix $C=9$ contacts per day and only adjust $a$. The results in Tables 6 and 7 also match up to the results in Table 3. Applying the agent-based approach to the intervention scenarios described in Section 3.2, we see that the results in Table 8 match with the results in Table 4 as expected, even maintaining the interesting behaviour around the timing of herd-immunity.

Now, using the agent-based approach, we look to explore and compare against the results presented in Section 3.5. It must be noted that in the deterministic setup we calculated the threshold of the maximum number of patients in hospital before an intervention is needed to guarantee that capacity is not breached. In this section, we want to calculate the number of Monte Carlo realisations that still result in a breached capacity using the calculated deterministic threshold and how the number of realisations that breach capacity decreases when that threshold is reduced. In particular, we look to do a more realistic investigation to the thresholds using the in-built stochasticity of the agent-based approach. This investigation can act as some sort of buffer for hospital management to understand how close to the threshold they can get before having to make changes. We note here that the agent-based approach tends to overestimate its deterministic counterpart when $\Delta t=1$, as approaches the deterministic solution from above when reducing $\Delta t$, as seen in Figure 16 in Appendix A. This implies that the results will be skewed towards a breach in the threshold. In this investigation we measure three metrics, the number of Monte Carlo realisations that go over $H_{max}$, the average maximum amount the realisations go above $H_{max}$ (to give an idea of how severely the threshold is breached), and the average maximum difference of the simulated hospital capacity $H$ from $H_{max}$ across all realisations, given a value of $H_{u}$ associated to a value of $\mathcal{R}_{0}$ from Figure 12. We measure these metrics against the percentage of $H_{u}$ which initiates an intervention, ranging from 90% of $H_{u}$ to 100% of $H_{u}$. We make several observations from the results in Table 9:

1. 1.

   The higher the value of $\mathcal{R}_{0}$, the higher the chance the system will be breached;

2. 2.

   The higher the value of $\mathcal{R}_{0}$, the lower the percentage of the threshold needs to be considered to reduce the chance of a breach;

3. 3.

   The higher the value of $\mathcal{R}_{0}$, the larger the system breach will be on average.

Intuitively, these observations make sense and should not necessarily come as a surprise. However, what is useful about this study is understanding what percentage of $H_{u}$ is needed to reduce the breaches. This result and process can be used by hospital administrations, planners and managers to prepare for an incoming wave and set aside surge capacity in case of a breaching realisation.