1. Introduction

The most recent outbreak overwhelming a healthcare system took place in west Africa in 2014-2015. A total of 28, 616 Ebola cases in Guinea, Liberia and Sierra Leone were reported to WHO as of June 10, 2016 with 11, 310 deaths [16], with cases occurring as far away as the United States.

An outbreak can overwhelm hospitals and clinics by affecting the community on a large scale. Healthcare facilities in developing countries are very often inadequate to fight against any widespread disease. Even in the United States, the healthcare system does not have enough infrastructure to fight against pandemics. According to the U.S. Department of Health and Human Services the total number of staffed beds in 2013 was 914, 513. Surprisingly, the number of staffed beds reduced about 38% from 1975 [23].

When a disease spreads rapidly, the resources required to treat patients effectively may fall short. However, resource limitations are often overlooked in modeling outbreaks. In the last few years, many researchers have used mathematical models to study Ebola virus disease (EVD), some of whom addressed the issue of limited hospital beds. In 2015, Drake et al. [9] investigated the 2014 Ebola outbreak in Liberia, where a dramatic increase in hospital capacity was observed during late August (to about 400 beds, an order of magnitude higher than a month earlier). They investigated scenarios in which hospital capacity was increased heterogeneously. They forecast 130, 000 cases by December 31 with the existing hospital capacity, but only 50, 000 for the same period when the number of beds was hypothetically ramped up to 1700. However, their model suggested that without rapid hospitalization, such an increase would not achieve containment. More recently, Njankou and Nyabadza (2018) studied the dynamics of EVD [15] using a modified SIR [deterministic] model with a time-dependent number of beds. With a limited number of beds, a backward bifurcation arose, complicating the control of EVD. Finally, Ahmad et al. (2016) used optimal control analysis on an SEIR model for EVD to study allocating resources among additional hospitalization, quarantine and vaccination components [2]. They compared constant vs. optimal control strategies in two cases: (ⅰ) hospitalization and quarantine and (ⅱ) hospitalization and vaccination, and concluded that optimal control is preferable for EVD. However, this requires that the healthcare system have enough resources for hospitalization and medication for a considerable number of susceptibles.

All these studies suggested increasing the number of hospital beds without bound, but each hospital has its own capacity limitations. When that limit is reached, hospital officials must decide how to respond to additional patients. This study aims to consider this decision in the context of the resulting disease burden on the community. We incorporate the idea of the healthcare system's carrying capacity (number of regular beds) $K$ to observe the way a disease spreads when the number of patients increases beyond the hospital's regular capacity. During an outbreak, if patients continue to come to the hospitals or clinics when those are already occupied fully, then hospital authorities have two options-they can either continue to admit patients (policy Ⅰ), or send the patients back as they lack the resources to care properly for them (policy Ⅱ). Here, we assume that in an overcrowding scenario a hospital can accept more patients than its carrying capacity, reducing infection in the community, but in return the quality of protection begins to break down, increasing both the death and infection rates within the hospital setting. If instead the hospital stops accepting more patients to preserve care quality, the rising infection rate in the community may lead to an epidemic. There is therefore a trade-off between the effects of high infection prevalence in the hospital vs. in the community.

Public health officials need to do some estimation while deciding their policy. They should adopt any policy to minimize loss. In an epidemic, however, estimating disease burden is not straightforward. During the epidemic, some infected people will die, and others will survive after a few days of suffering. It is very difficult to establish a relation between the loss related to deaths and loss due to the suffering of the survivors. Here we use the most commonly used idea of $DALY$s (Disability Affected Life Years) to calculate loss during an epidemic. This estimation has two components: $DALY$s for the survivals and $DALY$s for the non-survivals. This study uses dynamical systems models as tools to compare the consequences of both policies, in the aim of determining which policy by the hospital will better serve the community.

---

2. Mathematical model

In order to examine the effects of a hospital's admissions policy on disease spread, we develop a compartmental model with details focusing on within-hospital transmission. We begin by separating those in the hospital from those in the community, and then within each setting by infection status. While maintaining a simplified model, we pattern our outbreak on Ebola. Research shows nobody has been attacked a second time by the same strain of Ebola virus [12]. So, here we are considering an SIR model where we have two different susceptible classes- one is $S_c$, presenting healthy people living in the community, and the other is $S_H$, people admitted to the hospital prior to the outbreak, for reasons other than the epidemic disease. These last represent the individuals at higher risk of nosocomial infection when the hospital setting deteriorates in an overcrowding scenario. There are three different groups of infected people, $I_c$, $I_H$ and $M$, among which $M$ represents the people admitted to the hospital for different reasons and developed nosocomial infection with the disease of the epidemic before recovering from their primary diagnosis. The groups $I_c$ and $I_H$ (respectively) represent the people in the community with infection and people in the hospital with only the disease of our concern. Anyone from $M$ who recovers from the primary diagnosis moves to $I_H$, and those recovered from the secondary (epidemic) disease in $M$ move to $R_H$. At recovery, people move to $R_c$ from $I_c$ and to $R_H$ from $M$. Anyone from $R_H$ who recovers from the primary diagnosis moves to $R_c$. Here we consider $H = S_H+M+I_H+R_H$ and $C = S_c+I_c+R_c$.

In modeling infection rates, we use the standard incidence in the community as the community is saturated with people and so we introduce $\beta_c \frac{I_c}{C}$ as the per capita infection rate in the community. Now, the infected people in the community transfer (at rate $p$) to the hospital $(I_H)$ or move (at rate $\gamma_c$) to the recovered class $(R_c)$ or can die (at rate $d_c$) due to Ebola infection. For communities in which traditional burial practices introduce a significant source of additional infection, a term proportional to $d_c I_c \frac{S_c}{C}$ (the total death rate, multiplied by the proportion of the community available to be infected) is added to the base infection rate, but the net result is simply to increase the value of $\beta_c$. We assume that deaths occurring in the hospital are required to follow secure burial practices.

Since the population size is small inside the hospital, we consider mass-action incidence there. Hence, when the hospital is operating at normal capacity $(H\leq K)$, patients from $S_H$ either move to the nosocomial class at $\beta_H (M+I_H)$ or return to the community (at $q$) as susceptible. Patients from $M$ either transfer to $I_H$ (at rate $q$) or to $R_H$ (at rate $\gamma_H$) or can die (at rate $d_H$). Finally, patients move to $R_c$ from $I_H$ (at rate $\gamma_H$) and from $R_H$ (at rate $q$). Here, we assume the recovery rates from $M$ to $R_H$ and from $I_H$ to $R_c$ are equal, and the transfer rates (rates of transition between compartments) from $M$ to $I_H$ and from $S_H$ to $S_c$ are same. The death rates in the nosocomial class $(M)$ and in the infected class in the hospital $(I_H)$ are also assumed equal. We also assume that the epidemic will continue for a short term and consequently we ignore natural birth and death rates in our model.

With regard to the effects of hospital overcrowding, we assume a gradual deterioration in the quality of the hospital setting as admissions increase past the carrying capacity. To describe this deterioration, we introduce a function f defined by $f(H) = max(1, e^{\epsilon(H-K)/K})$ which amplifies the in-hospital infection rate and death rate (and slows down recovery) if the hospital continues to accept patients beyond its capacity $(H\geq K)$. While defining $f(H)$, we introduce $\epsilon$ as a calibration parameter to control the degree to which the situation deteriorates with overcrowding. The in-hospital infection and death rates increase by a factor of $f$ after admissions surpass the carrying capacity. Correspondingly, the in-hospital recovery rates ($q$ and $\gamma_H)$ are multiplied by $1/f$ when $H \geq K$. Finally, our model becomes:

$\begin{align*} \frac{dS_c}{dt} & = -\beta_c\frac{I_c}{C}S_c+\frac{1}{f}qS_H\\ \frac{dI_c}{dt} & = \beta_c\frac{I_c}{C}S_c-\gamma_cI_c-pI_c-d_cI_c\\ \frac{dS_H}{dt} & = -f\beta_H(M+I_H)S_H-\frac{1}{f}qS_H\\ \frac{dM}{dt} & = f\beta_H(M+I_H)S_H-\frac{1}{f}(\gamma_H+q)M-fd_HM\\ \frac{dI_H}{dt} & = \frac{1}{f}qM-\gamma_HI_H+pI_c-fd_HI_H\\ \frac{dR_H}{dt} & = \frac{1}{f}(\gamma_HM-qR_H)\\ \frac{dR_c}{dt} & = \gamma_cI_c+\frac{1}{f}(\gamma_HI_H+qR_H) \end{align*}$

where $d_c>d_H$

In our model, we ignore the case of people moving from $S_c$ to $S_H$ as during the epidemic relatively few people will be sick for reasons other than the epidemic, and during an epidemic people typically avoid going to hospitals [14]. Here, to simplify and focus analysis on in-hospital impact of overcrowding, we assume no hospital visitors are infected or contaminated, so no infection from hospital environment to community; and we assume that any contamination of healthcare workers leading to infection can be incorporated into the nosocomial infection term.

In our research, we have two additional differential equations to calculate the number of infections and number of deaths during the pandemic-

$\begin{align*} \frac{dI_c}{dt} & = \beta_c\frac{I_c}{C}S_c+f\beta_H(M+I_H)S_H \\ \frac{dD}{dt} & = d_cI_c+fd_H(M+I_H) \end{align*}$

where $I$ and $D$ represent the cumulative total number of infections and number of deaths respectively.

Now we define our cost function which is determined by the number of infections and by the number of deaths and will estimate the total burden of the disease. The accounting breakdown of this cost function is shown in Fig. 2. It can also be considered as a burden function or loss function. This cost function will determine the number of $DALY$s which can be thought of as loss of healthy life. WHO formulate how to calculate this loss [17]. As WHO prescribed, the $DALY$ is composed of Years Lost due to Disability $(YLD)$ resulting from infections and Years of Life Lost $(YLL)$ caused by premature death. To estimate those values, we have $YLD = I\times DW\times L_1$ and $YLL = D\times L_2$, where $L_1 =$ average duration (in years) of the case until remission or death (the duration of illness), which is the reciprocal of the recovery rate (for survival) or of the death rate (for non-survival) and $L_2 =$ standard life expectancy at age of death (in years) = average life expectancy of the community － average age at premature death. The term DW is the disability weight for the disease which ranges from 0 (perfect health) to 1(death). It can be thought of as the proportional reduction in perfect health due to any adverse health condition. So, the total burden of disease can be represented as the sum of $YLD$ and $YLL$, which leads us to $-$

${\rm{Total}}\;\;DALY = I \times DW \times {L_1} + D \times {L_2}$

But, $L_1$ will be different in the community and in the hospital as well as different for the survival and non-survival. So, to make our analysis more accurate we divide $I$ in four subgroups (Fig. 2).

Then, finally our cost function is-

$J(t) = (I_{SN}\times \frac{1}{\gamma_c} I_{SH}\times \frac{1}{\gamma_H} I_{NH}\times \frac{1}{d_H} I_{NN}\times \frac{1}{d_c}) \times DW +D\times L_2$

where the numbers of people in ${I_{SN}},{I_{SH}},{I_{NH}}$ and ${I_{NN}}$ classes are tabulated by integrating the respective exit rates $( \gamma_c I_c, \gamma_H (I_H+N), d_H (I_H+N)$, and $d_c I_c)$.

In many disease outbreaks, patients who recover and are discharged from the hospital may return shortly afterward from complications caused by the illness. We therefore also considered an extension of the model described above, in which individuals in $R_c$ move either to $R_H$ (readmission) or to a permanently recovered class $R$. However, a readmission rate of 15$%$ made an insignificant change (less than 0.04$%$) in the final outbreak size reported in the analysis section. So, in the remainder of the study we ignore the impact of readmissions during the outbreak.

---

3. Parameter estimation

In our work, we assume the population in the community is 100, 000; among those 99, 980 are susceptible and 10 are infected. We also assume at the beginning of epidemic 10 people were already in the hospital for different reasons other than the epidemic. We know the number of hospital beds per 1000 people in Sierra Leone is 0.4 [1]. Using this documented data, we assume a hospital in the area has a capacity of 40 beds. Here we ignore the number of health workers while calculating the number of beds in hospital.

While estimating parameters, we try to take the values from the same epidemiological context (Sierra Leone, 2014) to make our analysis more appropriate. The infection rate in the community $( \beta_c)$ is taken as 0.455/day [4]. Here, we take the reciprocal of the mean of boundary values of the range. Unfortunately, we did not find any documented data for the infection rate in hospitals. However, in our literature review, we find that control measures (including education and contact tracing followed by quarantine) reduce the infection rate by 50$%$ (from 0.38/day to 0.19/day) for the 1995 Ebola epidemic in Congo and by 73$%$ (from 0.33/day to 0.09/day) for the 2000 Ebola epidemic in Uganda [7]. So, we assume the transmission rate in our analysis is reduced by 61.50$%$ in the hospital (mean of 50$%$ and 73$%$). To get infection rate in a hospital, we reduce our value of $\beta_c$ by 61.50$%$ and then divide it by the carrying capacity as we consider mass-action incidence inside the hospital. Here, we divide by the carrying capacity assuming a hospital will be fully occupied with patients after very few days of the emergence of the disease. Finally, we get 0.004375/person-day as infection rate in hospitals. We also have 0.057/day as recovery rate in hospitals $( \gamma_H)$ [5]. Here, we use the reciprocal of the sum of days for onset of symptoms to hospitalization and days for hospitalization to death. However, we did not find any well documented data for the recovery rate in the community. So, here we use the data from a research work on Ebola virus infection in rhesus monkeys to compare recovery rates with and without treatment. We take the ratio of the survival periods of untreated monkeys (8.3 days) to treated (11.7 days) and apply it to the recovery rate in hospitals to estimate the recovery rate in the community ($\gamma_c = (8.3/11.7)\times \gamma_H = (8.3/11.7)\times 0.057/$day = 0.04/day) [10]. We also did not find any authentic source for the average length of stay in hospitals in Sierra Leone and not even in Liberia or in Guinea. Here, we consider the average length of stay in hospital for Uganda and take the reciprocal to calculate the value of parameter q which is 0.067/day [11]. For the transfer rate of infected people from community to hospitals $(p)$ we take the weighted mean 0.184/day of some earlier documented estimation (Table 1).

We use the duration (mean of two values: 5.0 and 6.6) from onset of symptoms to death and take the reciprocal to estimate death rates. Thus, we have 0.172/day [3] and 0.102/day [22] as death rates in the community and in hospital. These two values are taken as the reciprocal of the number of days.

A careful literature review turned up almost no data on how overcrowding deteriorates the quality of the hospital setting. During the 2014 Ebola outbreak, the WHO reported on one clinic in Liberia with 120 beds which had admitted as many as 210 patients, 75% more than its carrying capacity [18], and noted that when a new 20-bed clinic opened in Liberia's capital city, it was immediately overwhelmed with more than 70 patients [19], although the report does not indicate how many patients were able to remain there. Likewise, a MSF team rehabilitated a 40-bed facility elsewhere in Liberia to the point that, two weeks later, it was caring for 137 suspected Ebola patients [14], although the report does not specify the nature of the rehabilitations, which may have included expansion or permanent new beds. In line with the first, more moderate figure above, we consider in our model the case where a hospital can admit up to a maximum of 70% more (by doing arrangements on floor and establishing temporary tents inside the hospital premises) than its carrying capacity [13]. In the absence of further data, we assume that at maximum overcrowding the quality of the hospital setting degrades by 50%, i.e., the infection and death rates in hospital increase by 50%. This assumption leads to a value of $\epsilon =$0.48067 in defining the function $f(H)$. Although the numerical results below clearly follow from this estimate, using different values for $\epsilon$ yields qualitatively similar results.

For the estimation of the value of cost function we have the average age of infection in Sierra Leone as 28 years [20] and the average life expectancy in Sierra Leone as 57.39 years [6] which gives $L_2 =$ 57.39 - 28 = 29.39. Unfortunately, there is no documented data for the $DW$ for Ebola. However, we find a very well documented study on Global Disease Burden where $DW$ values are reported for three types of infectious diseases- mild, moderate and severe. For Ebola, we choose the severe case and the $DW$ for that is 0.133 [21].

---

4. Analysis

To find the equilibria of our dynamical system we set all the equations of our model equal to zero and solve those. After doing some calculation, we get one solution set $(I_c = S_H = M = I_H = 0, S_c \neq 0, R_c \neq 0, S_c+R_c \leq C(0)+H(0))$ which represents an infinite number of disease-free equilibrium points. Because these equilibria are non-isolated (in two dimensions, in fact), their stability cannot be analyzed. In general, solutions of the system will eventually approach one of them, first rising to an epidemic peak if the control reproduction number¹ $R_0>1$. Further qualitative analysis of the system is therefore limited to deriving $R_0$, which remains crucial to understanding the behavior of the outbreak. So, we use the Next Generation Method [8] on our model to derive the expression for the control reproduction number, $R_0 = \frac{\beta_c}{\gamma_c+p+d_c}$, which gives 1.149 for our baseline estimation of control reproduction number and for the basic reproduction number (setting the transfer rate p to 0) it gives 2.146. In the analysis that follows, we use $R_0$ as an epidemiological index by which to compare the effects of the two admission policies.

¹ The control reproduction number is distinguished from the basic reproduction number by the inclusion of control methods. Here we use the familiar notation $R_0$ to refer to the control reproduction number.

In our numerical analysis, we found the hospital will reach its carrying capacity in 26 days and the maximum 1.7 times its carrying capacity in 37 days (Fig. 3a). We also found the epidemic will continue for 136 days and 128 days when policy Ⅰ and policy Ⅱ are adopted respectively. In estimating the cost function, we found a loss due to the epidemic of approximately 2.35 million $DALY$s in either case, but policy Ⅱ (no hospital overcrowding) results in over 500 fewer $DALY$s than policy Ⅰ.

Although Policy Ⅰ results in fewer infections, it leads to more deaths and thus a higher overall disease burden (Table 3). Thus, the policy of admitting patients only up to a maximum of the hospital's carrying capacity is better (by 507 $DALY$s) for our set of parameter values. Stopping admissions after reaching the carrying capacity thus produces a lower disease burden in this case.

However, the epidemic will be longer if the hospital works under policy Ⅰ (Fig. 3b). This happens because health care benefits in the hospital slow down the epidemic. Our set of parameter values shows how the hospital's admission policy in terms of admitting patients affects the burden of the epidemic at a baseline level.

To find which parameters most influence the final cost of the epidemic, we perform a sensitivity analysis on our model. In this process, we increase the value of all the parameters by 1$%$ to see the effect of these changes on the final value of the cost function. Our analysis (see Fig. 4) shows infection rate ($\beta_c$), recovery rate ($\gamma_c$) and the death rate ($d_c$) in the community impact the total loss an order of magnitude more strongly than all the other parameters. So, we vary the infection rate, death rate and recovery rate to observe the behavior of our system.

To check the effect of infection rate on the impact of the two policies, we vary the infection rate and try to establish a relation between the infection rate and the entire loss (value of the cost function) due to the epidemic. Here, we used the range 0.380/day to 0.515/day for the value of the infection rate in the community ($\beta_c$). We assume the infection rate in the hospital ($\beta_H$) will change at the same ratio as the infection rate in the community. In Fig. 5a, the difference in behavior of the policies is manifest near the $R_0 = 1$ (indicated by the up arrow) bifurcation though away from the bifurcation it is obscured due to the issue of scale. Then, we generate another graph (Fig. 5b) showing the difference in total burden ($DALY$s) associated with the two policies with the same varying infection rate in the community ($\beta_c$). It clearly manifests the behavior difference due to the policies everywhere. In Fig. 5b, the $DALY$ difference is taken as the difference between $DALY$s associated with policy Ⅰ and $DALY$s associated with policy Ⅱ. So, policy Ⅰ is better when $\beta_c$ ranges from 0.401/day to 0.438/day ($R_0\gtrsim$1, here $1.012\leq R_0\leq 1.106)$ and policy Ⅱ is better when $\beta_c$ is higher than 0.438/day ($R_0\gg$1). Here, the range for $R_0$ does not begin just after 1.00 as the epidemic is not big enough to fill up the hospital unless the value of $R_0$ reaches 1.012. However, the impact of the epidemic is the same regardless of the policies if $\beta_c$ is 0.396/day or less $(R_0\leq 1.012)$, since the hospital never fills up.

Then, we try to observe the change in the behavior of the burden of the epidemic if the death rate is changed. Here, we also assume the death rate in the hospital $(d_H)$ changes at the same ratio as the death rate in the community $(d_c)$ is changed. To see the effect of change in death rate we use the range 0.070/day to 0.320/day for the value of $d_c$. In Fig. 6a, we have the total burden of the epidemic with varying death rates. However, the behavior difference away from the $R_0 = 1$ (indicated by the up arrow) bifurcation and near the horizontal axis is not manifest. Then, we generate a graph (Fig. 6b) showing the difference in total burden, where values below the horizontal axis imply policy Ⅱ generates more $DALY$s and values above the horizontal axis imply policy Ⅰ generates more $DALY$s. It is evident from Fig. 6b that policy Ⅰ is better when the death rate (in the community) ranges from 0.180/day to 0.231/day $(R_0\gtrsim 1, 1.013\leq R_0\leq 1.126)$ and policy Ⅱ is better when death rate (in the community) is 0.180/day or less $(R_0\gg1)$. However, both policies result the same when the death rate (in the community) is 0.231/day or above $(R_0\leq 1.013)$.

The same pattern is found while varying the recovery rate $\gamma_c$ from 0.028/day to 0.11/day(Fig. 7). Policy I is better when the recovery rate (in the community) ranges from 0.05/day to 0.09/day $(R_0\gtrsim 1, 1.02\leq R_0≤ 1.12)$ and policy Ⅱ is better when this rate is 0.05/day or less $(R_0\gg 1)$. However, no impact on the behavior of the policies is observed when the recovery rate (in the community) is 0.099/day or above $(R_0\leq 1.02)$.

Our simulations produce an overall fatality rate of about 80$%$, based on parameter estimates from the Ebola epidemic in West Africa in 2014, which is significantly different from the actual rate in that epidemic (about 50$%$). This happens because we used published data (which included overcrowding) to derive baseline (non-overcrowding) estimates and then applied overcrowding on top of that. Simulations using other parameter values produced qualitatively similar results.

Our results are dependent on the choice of scaling parameter $\epsilon$ and maximum possible overcrowding inside the hospital. Changes in these two values will either widen or narrow the breadth of the interval of $R_0$ for which overcrowding a hospital is better. However, the qualitative results will remain the same although the upper $R_0$ threshold and the magnitude of the $DALY$ difference will vary.

---

5. Discussion

It is not surprising to imagine that the effect of any epidemic will be worse if hospitals stop admitting patients after reaching carrying capacity. However, our work shows some interesting results. Our simulations of an Ebola epidemic indicate that it is sometimes better to stop admitting patients after hospitals reach their carrying capacity, rather than continue to admit patients and overcrowd the facility (in the latter case we assumed a hospital can choose to continue to accept patients up to a maximum of 70% more than its carrying capacity). There is a narrow window where $R_0$ is slightly greater than 1, where overloading a hospital is better because this is enough to handle most cases. When the system is overwhelmed anyway ($R_0\gg 1)$, overcrowded hospitals spread infections and, in that scenario, maintaining hospitals at their carrying capacity is better, because the quality of the hospital setting becomes so compromised when highly overcrowded that it does more to incubate infections than to limit them. However, when $R_0\leq1.012$, a hospital's admission policy beyond its carrying capacity does not matter because there is no epidemic large enough to fill the hospital completely.

The relation between the policy of admitting patients beyond hospitals' carrying capacity and burden of Ebola is not straightforward: it depends on the value of $R_0$. It is better to continue admitting patients beyond their carrying capacity when $R_0 \gtrsim$1 and better to stop admitting after reaching carrying capacity when $R_0$$\gg$1 in order to avoid heavily overcrowded hospitals increasing infection risk. However, the $R_0$ value delineating these two cases depends on how many more patients hospitals are able to admit. Different choices of scaling parameter $\epsilon$ and maximum possible overcrowding inside the hospital will change the interval of $R_0$ values on which the choice of policy depends, but not the qualitative result that slightly overwhelmed healthcare facilities can continue to protect the community, while grossly overloaded ones risk becoming sources of infection as crucial protective measures break down.

Our decision to ignore (or suppose negligible relative to other sources) infections transmitted in-hospital to visitors potentially underestimates the total epidemic size but allows analysis to focus on the two distinct epidemiological settings.

Future work can be done to investigate whether or not this result is true for any infectious disease epidemic. Here, we deal with a non-vector-borne disease. It will be interesting to see how the two policies behave when a vector-borne disease is taken in account. Further studies could also investigate the impact of other resource limitations such as antiviral stockpiles.

---