Modelling strategies to organize healthcare workforce during pandemics: Application to COVID-19

Protection of the healthcare workforce is of paramount importance for the care of patients in the setting of a pandemic such as coronavirus disease 2019 (COVID-19). Healthcare workers are at increased risk of becoming infected. The ideal organisational strategy to protect the workforce in a situation in which social distancing cannot be maintained remains to be determined. In this study, we have mathematically modelled strategies for the employment of the hospital workforce with the goal of simulating the health and productivity of the workers. The models were designed to determine if desynchronization of medical teams by dichotomizing the workers may protect the workforce. Our studies model workforce productivity and the efficiency of home office applied to the case of COVID-19. The results reveal that a desynchronization strategy in which two medical teams work alternating for 7 days increases the available workforce.


Introduction
During the spread of worldwide pandemics, protecting and supporting caregivers to maintain the workforce in hospitals is a crucial and challenging task. Under such extraordinary situations, the efficiency of the healthcare workforce is threatened by several factors including: 1) infected patients, 2) infected co-workers and 3) infected persons outside the hospital. This is of high importance in the current coronavirus disease 2019 pandemic  given that the virus may be transmitted by asymptomatic persons and the long incubation period averaging 5 days (Guan et al., 2020;Zhou et al., 2020;Lauer et al., 2020). Unlike for most other professions, social distancing is typically not possible in medical teams where healthcare workers are required to work in close contact with patients and coworkers. Currently, little is known about which work organizational strategies in the hospital are most suitable for protecting the healthcare workforce. Potential strategies may include prolonging duty hours and thereby limiting interactions between coworkers or complete desynchronization of the workforce in which teams are dichotomized and each half of the team work alternate weeks.
In this report, we perform simulations of biophysical models of coronavirus epidemics of COVID-19 in a healthcare working team, and discuss the efficiency of different work strategies during the viral outbreak.
Different types of mathematical models have been used to study epidemiology. In the classical Susceptible-Infectious-Recov ered (SIR) model, a susceptible worker (S) can be infected (I), and the infected person can recover (R), without the risk of reinfection. Variants of this model are the Susceptible-Infectious-Susceptible (SIS) and Susceptible-Exposed-Infectious-Recovered (SEIR) models (Getz et al., 2018) and the Susceptible-Infected-Recovered-Deceas ed (SIRD) models (Sinha et al., 2020). The SEIR model in which the recovered patients are susceptible again, has already been used for modeling purposes during the COVID-19 outbreak (Rocklöv et al., 2020). Also, SEIR models have been used to model control of expansion in the context of the COVID-19 and have been extended to include patients age and asymptomatic cases (Pan et al., 2019;Prem et al., 2020). Furthermore, SEIR models have been applied in this context by including persons in quarantine (QSEIR model) Peak et al., 2020). A model based on Microscopic Markov Chain approach have been proposed to study the effect of confinement (Arenas et al., 2020).
Further extensions have recently been proposed to better model COVID-19 transmission dynamics. Mobility data from the United States have been used to create social-distance metrics and combined with epidemiological models to assess infection rates (Badr https et al., 2020). Recent work has analyzed the impact of population heterogeneity in epidemic waves (Neipel et al., 2020). Models of the transmission dynamics in Switzerland have been developed to study the effect of the implementation of non-pharmaceutical interventions (NPIs) on time-varying reproduction number (Joseph et al., 2020) as well as to explore counterfactual scenarios of earlier and later implementation of NPIs (Althaus et al., 2020). Optimization of resource sharing of multiple hospitals has been studied to optimize the use of ICU beds (Lacasa et al., 2020). However, the question on how to organize healthcare units during pandemics remains unclear.
Here, we put forward SIS and SIR models by dividing the infectious persons into a latent and infected state in order to account for the potentially long asymptomatic phase of COVID-19. Next, have developed a family of time-dependent compartmental models with and without reinfection by adapting the SAIR model (Althaus et al., 2010;Leung et al., 2018;Kahn et al., 2019). We added a variable to account for work W and built mathematical models for COVID-19 described by ordinary differential equations (ODE). The family of models is denoted Susceptible-Latent-Infec ted-Recovered-Susceptible-Work (SLIRSW) model. In our SLIRSW models, the workers can be healthy and susceptible to infection, S, infected in the incubation period, L, infected presenting symptoms, I, and after recovery temporarily immune to new infections, R. Eventually the recovered workers lose their immunity and become susceptible again. Additionally, the models are adapted to account for work, W, done by the non-infected workers. There are two limiting cases of particular interest: when the recovered do not gain immunity (SLIW) and when the immunity is permanent (SLIRW).
To investigate the possible scenarios for workforce organization, we use ODEs describing the dynamics of the models including time-dependent parameters in which the rate changes based on their location, that is, in-hospital compared with home office. Deterministic models may not be accurate to describe the dynamics when small populations are present (McKane et al., 2005). In the context of COVID-19, stochastic epidemiology models have been used to study confinement strategies (Bittihn et al., 2020).
For this reason, we studied the limitations of the deterministic model by reformulating them in terms of a master equation (Gillespie, 1977), and we determined the range of parameters where stochastic models were necessary. Later, the stochastic models are used to investigate space separation of the workers.
The rest of the report is organized as follows: we formulate mathematical models to investigate organizational strategies of the hospital workforce and their effect on productivity during a pandemic (Section 2, and present the results of simulations (Section 3). Finally, our results are summarized and discussed (Section 4).

Models
In this section we present the mathematical models used to explore the different strategies to organize healthcare workers.
Our goal is to create a model to study the pandemic from its beginning, where the number of infected workers is 0, to later time points where the population becomes immune, or when the virus becomes endemic. Therefore we are studying both, the steady state and the transient dynamics.

Sketch of the models
In the normal/control strategy, the SLIRSW model describes the transmission dynamics among one group of healthcare workers that follows their usual schedule at the hospital. In that scenario, workers can be infected by coworkers and by the general population. This model is described in detail in Section 2.2.1. This normal strategy is compared with a desynchronized model, in which the healthcare workers are divided into two non-interacting teams of equal size. At any time, one team is at the workplace and the other team stays at home. A dichotomous work switch is implemented as follows ( Fig. 1A): during the first week Team 1 stays at the workplace and Team 2 at home; during the second week Team 1 stays at home and Team 2 at the workplace; next the teams continue alternating their location at the end of each week. When at the workplace, each healthy susceptible worker can become latent at an infection rate that is larger than when staying at home. Latent workers develop the infection at the workplace, at the same rate as at home. Infection of a susceptible by a latent worker can only occur in the workplace. We also assume that infected workers recover at the same rate at the workplace than at home, becoming recovered and immune and eventually losing their immunity. This model is described in detail in Section 2.2.1. Fig. 1 shows the SLIRSW model with compartmentalization.

Mathematical formulation
In this section we present the mathematical details of the SLIRSW model.

Normal strategy: One group of healthcare workers
In the normal strategy, our infection model describes the evolution in time t of the number S of healthy workers susceptible to be infected, the number, L, of infected workers in a latent state, the number, I, of infected workers presenting symptoms and therefore isolated, and the number R of recovered workers who gained temporal immunity and can not be infected. The dynamics of each population are described by the following system of ODEs: where a is the infection rate coming from the general populatione.g. infected persons in the city where the hospital is located-, resulting in an infected worker in a latent state, b is the infection rate when the person is infected by a co-worker, s is the inverse of the incubation period, or the activation rate, at which a patient in a latent state starts presenting symptoms. r is the recovery rate, where an infected patient recovers and becomes healthy, and m is the loss of immunity rate at which a recovered patient becomes susceptible.
By considering a to be constant, we assume that there is a roughly constant fraction of infected individuals among the general population that becomes infected and can infect the healthcare workers. Although more complex infection rate functions that introduce an infection curve could be considered, we use a constant force of infection for simplicity to represent infections from the general population. One interpretation of this approximation is to consider it as a worst-case approximation when the infection rate a is fixed at its highest during the whole epidemic. This worst-case approximation would produce a lower bound on the number of healthy workers.
W is a function that keeps track of the work output of the workers, for simplicity we assume its growth is proportional to the number of available workers. For simplicity, we consider the death rate negligible and we do not consider recoveries coming from latent patients. In our simulations reported below, we set the initial condition to Sð0Þ ¼ 300; Lð0Þ ¼ 0; Ið0Þ ¼ 0; Rð0Þ ¼ 0, and Wð0Þ ¼ 0, thereby representing a large unit with 300 workers.
As a function of m, there are two limiting cases of this model. When m ¼ 0, the immunity gained is permanent, we denote this model as SLIRW. When m ! 1, the time to go from R to S tends to 0, we will refer to this model as SLIW.

Desynchronization strategy: Two groups of healthcare workers
Our main goal is to identify strategies to protect healthcare workers. The first organizational strategy we present is a desynchronization strategy that we compare with the normal strategy presented in Section 2.2.1. The desynchronization strategy consists of the division of the workers into two groups, one group works for one week in the hospital while the other works at home, and the week after the two groups swap their roles. Accordingly, we model it as two groups of healthy persons susceptible to be infected, S 1 and S 2 , two groups of infected persons in a latent state, L 1 and L 2 , two groups of infected persons presenting symptoms, I 1 and I 2 , and two groups of recovered R 1 and R 2 . The dynamics of each population are described by the following system of ODEs.
for the first group and dS 2 dt for the second group. The rate of change of the work output obeys For this model, the parameters that are purely related to the disease, r; s and m remain the same, while the parameters that determine the epidemiological spreading vary. Accordingly, we simulate our model for the first seven days with b 2 ¼ 0 because the workers are not in contact with the other co-workers when they do home office. Additionally we define a H as the value a 1 ðtÞ (a 2 ðtÞ) takes during the home office of group 1 (group 2). We consider a H < a because we assume that the probability of becoming infected by the general population is higher when working at the hospital. Every seven days, we switched the group's roles, by exchanging their parameter values. We assume that the productivity decreases during home office, therefore we choose AðtÞ ¼ 1 and BðtÞ ¼ W H during the odd weeks and BðtÞ ¼ 1 and AðtÞ ¼ W H during the even weeks. Then, a i ðtÞ; b i ðtÞ; AðtÞ and BðtÞ can be written as Our model assumes that b remains constant when the desynchronization strategy is applied. The rational for this assumption is that with half of the workers in the same space, the time for them to interact would double. However, some procedures (e.g. surgical interventions) would require the workers to interact, increasing the value of b. The current lack of data does not allow us to estimate b during desynchronization, and the assessment of the true value is yet an open question.

Spatial divisions
Another strategy we consider, in addition to desynchronization, is the possibility of dividing the space in the hospital into spatial compartments.
In order to model space division, we divide our variables S i ; L i ; I i ; R i , in non-interacting compartments, S ði;jÞ ; L ði;jÞ ; I ði;jÞ ; R ði;jÞ with i ¼ 1; 2; j ¼ 1; . . . ; N and N is the number of compartments. Since the space is divided by N where N is the number of compartments, we assume the rate b is multiplied by N, that is, b ! Nb.

Stochastic model
In this section, we present a stochastic formulation of our model in terms of a Master equation, which we use to study smaller units, e.g. 75 workers instead of 300. Here, due to the small size, stochastic effects can become dominant.
We consider the susceptible, S, latent, L, infected showing symptoms I, and recovered R, coupled with the following reactions S can be infected by a non-co-worker resulting into a L with a rate T 1 . See Table 1. S can be infected by a co-workers resulting into a L with a rate T 2 . See Table 1. A latent infected worker L can start presenting symptoms and be sent to quarantine at rate T 3 . See Table 1. An infected worker I can recover and become R at rate T 4 . See Table 1. A recovered worker R can lose their immunity and become S at rate T 5 . See Table 1.
The SLIW and SLIRW models are reformulated in terms of the following Master Equation @PðX; tÞ @t ¼ X j ðT j ðX À r j ; tÞPðX À r j ; tÞÞ À T j ðX; tÞPðX; tÞÞ Where PðX; tÞ is the probability of being at state X at time t, where X ¼ ðS; L; I; RÞ in the control model and X ¼ ðS 1 ; L 1 ; I 1 ; R 1 ; S 2 ; L 2 ; I 2 ; R 2 Þ in the desynchronized model. The transition rates, T j , and the state changes, r j , are described in Table 1.
The Work done, W, is described as the ODE of the deterministic models, Eq. (14).

Desynchronization increases the available workforce during the peak of infected healthcare workers
In this section we performed numerical simulations of our dynamical systems model to study the effect of the desynchronization. Our strategy consists of comparing simulations of the model with desynchronization and the model without desynchronization in terms of number of healthy workers, and amount of work done.

Simulations of the SLIW model
First, we study the effect of the infection rates on the health of the workers. In Fig. 2 we show the simulations of our systems of ODEs with two representative sets of parameters, low infection rate (Fig. 2Ai) and high infection rate (Fig. 2Aii). Regardless of the values, we observe that the strategy of dividing the workers into two groups decreases the number of infected workers, thus increasing the available workforce. In Fig. 2B we show the number of workers as a function of the infection rates. We observe that with the desynchronization strategy the number of healthy workers is higher than with the normal strategy. So far, we have observed that the desynchronization strategy increases the number of healthy workers. Because a desynchronization strategy includes decreased productivity of home office, we next determined the overall economic impact of this strategy. Therefore, we simulate our model for 6 weeks and we represent the value of W divided by the total number of workers as a function of a and b, for different values of home office efficiency (Fig. 2C). We observe that if the productivity decreases by 50% during the home office week, the normal strategy performs better than the desynchronization strategy. For higher values (75% or 100%) of the home office productivity, the desynchronization strategy outperforms the normal strategy for high infection rates.

Simulations of the SLIRW model
Next, we study the effect of the infection rates in a model that does not allow for reinfection. In Fig. 3A we show the simulations of our systems of ODE with two representative sets of parameters. Regardless of the values, we observe that the strategy of dividing the workers into two groups increases the number of healthy workers when the maximum of infected healthcare workers is reached. In Fig. 3B we show the number of minimum healthy workers as a function of the infection rates during a 6-week simulation. So far, we have observed that the desynchronization strategy also increases the number of healthy workers in the situation without reinfection. To determine the economic impact, we simulate our model for 6 weeks and we represent the value of W as a function a and b, for different values of home office efficiency (Fig. 3C). The desynchronization strategy outperforms the normal strategy for high values of home office efficiency. However, to achieve the same amount of work done, the desynchronization strategy requires higher values of home office productivity than in the SLIW model.

Simulations of the SLIRSW models
In this section, we study the remaining scenario, where the worker gains immunity after recovery, which is however not permanent. We consider that this immunity is lost at a rate m which can vary from 0 (SLIRW model) to 1 (SLIW model).
In Fig. 4A we illustrate a 12 week simulation of the number of healthy workers for different values of m. They showed similar dynamics before the peak of infected healthcare workers, all of them reaching a similar minimum value. Thereafter, each condition tends to a different steady state, which decreases with the value of m.
In Fig. 4B we show how the minimum value and the steady state vary as a function of the parameter m. We do not observe remarkable differences in the minimum value, and desynchronization reduces the number of maximum infected workers for all values of m. The steady state increases for values close to 0. Therefore, if the immunity is preserved for a long time, the system tends to the SLIRW results.
In Fig. 4C we show the WðtÞ for multiple values of m. During the first weeks we do not observe differences. After that, the lower the value of m, the higher the value of WðtÞ. This effect is due to the fraction of non-healthy workers after reaching the peak of infected healthcare workers (Fig. 4A). In this scenario, we observe that the desynchronization strategy is beneficial for high values of m and only detrimental for low values of m for long periods of observation.
In summary, our results show that the effect of desynchronization can be helpful to protect healthcare workers for all values of m, and the dynamics for all values of m are in between the extreme cases SLIW and SLIRW. 3.1.4. Impact of home productivity on the performance of the two strategies for the limiting cases Next, we compared the SLIW and SLIRW. To this end, we defined i R as the infection rate and we computed the value of W as a function of the home office efficiency for different values of the infection rate, by choosing a ¼ 0:1i R and b ¼ 0:001i R (Fig. 5).
We observe that the SLIW model requires a lower home office efficiency than the SLIRW model (Fig. 5A,B). Next, we computed the minimal value of home office efficiency needed to not have loss in work using the desynchronization strategy (Fig. 5C).
Note that we have fixed a ratio of a=b to compare the SLIW and the SLIRW model, and a different value could change the quantitative behavior of what we observed in Fig. 5. However, this selection is a representative value of the qualitative dynamics which allows for the understanding of the difference between the extreme cases.

For small values of a and high values b the stochasticity is dominant
Our deterministic model allowed us to identify differences between the normal and the desynchronization strategy. However, it remains unclear if some features could be hidden because of stochasticity. In this section, we study the stochastic effects on our model by performing Gillespie simulations.
In Fig. 6A, we show that for low levels of a and low values of b, the stochasticity does not play an important role due to the low number of infections. For high values of a, regardless the value of b, stochasticity does not play an important role, in that situation, stochastic simulations are oscillating around the mean-field dynamics (6C, D). However, we found that for low rates of a and high rates of b, the stochastic effects are dominant (6B). This result led us to study them further as shown below.

Spatial desynchronization introduces a purely stochastic protective mechanism
In this section we explore the effect of dividing the work space by reducing the number of persons working in a team. The workers are assigned to one small team and are never in contact with workers from other teams. In doing so, we simulate the stochastic model with spatial compartmentalization described in Section 2.3.
In Fig. 7A,B we observe that several stochastic trajectories remain unchanged (i.e. at 1), indicating that all individuals stay non-infected in those simulations. The number of non-infected trajectories increases with the number of compartments, while the mean field dynamics remain the same. In Fig. 7C we observe that, on average, an increase in the number of compartments is translated into a higher number of healthy workers. Note that the changes in the number of compartments come with an increase in b i resulting in the ODE system giving the same result (Fig. 7A). Therefore, the spatial division incorporates a purely stochastic protective mechanism to our system. This protective mechanism points to possible problems with the interpretation of certain regions of Fig. 2B,C and Fig. 3B,C.
Specifically, the upper left region of the heatmaps (high b, low a).
In that regime, the mean-field dynamics possibly overestimate the number of infected workers.

Discussion
Given the urgent need to protect caretakers, we present here modelling approaches that address the impact of desynchronization of healthcare workforces. We have developed a family of time-dependent compartmental models (SLIRSW) with and without reinfection by adapting the SAIR model of Althaus et al. (2010), Leung et al. (2018) and Kahn et al. (2019). A timedependent compartmental model was included to account for the desynchronization of healthcare teams. In addition to the availability of healthcare workers, we modelled productivity by incorporating different levels of work performed at home office.
First, we studied the limiting case where workers turn susceptible again immediately after recovery (SLIW model), which showed that the desynchronization strategy is associated with an increase in the number of healthy workers compared to no desynchronization. This effect is present with both high and low levels of infection rates (Fig. 2).
Next, we incorporated productivity of the workers in the period of home office. In practice, productivity of home office may depend on the tasks that may be done outside the hospital such as writing reports. We have determined the overall productivity as a function of infection rates and home office efficiency. Figs. 2 and 5 show that, for our case study model of COVID-19, a decrease of productivity to 50% during home office (only half of the time can be used for productive work) does not imply a substantial decrease of the overall productivity. However, if the productivity at home office is above 75%, overall productivity is increased with a desynchronization strategy for high infection rates.
Second, we studied the case, in which the workers gain permanent immunity to the disease after recovering (SLIRW model). In this SLIRW model, the protective effect of desynchronization reduces over time (Fig. 3). In this situation, however, the number of healthy workers with desynchronization increases, especially during the peak of infected healthcare workers (Fig. 3). This observation is potentially the consequence of workers being immune to the infection thereby not requiring a desynchronization anymore. The protection of the healthcare workforce (Fig. 3B) and productivity (Fig. 3C) also depends on the infection rate in the SLIRW model.
Third, we studied an intermediate scenario in which the worker gains immunity after recovery, which can be lost after some time. We quantified the effect of the protective measure (Fig. 4), which is between the two extreme cases of SLIW and SLIRW.
Next, we aimed to determine the ideal level of productivity for each model (Fig. 5). The ideal home office productivity strongly depends on infection rates in the model with reinfection (SLIW). In this model productivity is required to be around 60% for infec-tion rates 3 and 10. In the model without reinfection (SLIRW), which is rather the case with COVID-19, the home office productivity needs to be higher in order to keep overall productivity similar between the one and two group strategy.
However, the way the model was designed, we artificially included a handicap to the desynchronization strategy. As long as it remains unclear if reinfection is possible, the desynchronization strategy would be maintained. Such strategy would be stopped for recovered workers if reinfection does not occur because of immunity. A possible future direction would be the study of optimization strategies to manage healthcare workers with adaptive strategies that analyze the temporal immunity of the recovered workers. Next, we reformulated the model in terms of a master equation to study the effect of the intrinsic noise. In Fig. 6 we showed that for low levels of a and high levels of b, the stochastic effects can be dominant, and therefore, in that regime a deterministic model would not be an accurate approximation of the process.
The stochasticity was of special importance in our next in silico experiment, where we studied spatial separation, due to the low number of individuals in each compartment. We have shown ( Fig. 7) that for low levels of a i , the spatial separation reduces the spread of the disease among the workers. We investigated the origin of this effect by studying the distribution of the number of infected workers, and we have shown that most of the simulations had no infections (Fig. 7A,B), reducing considerably the average number of infected workers (Fig. 7C). This stochastic protective measure is a discrete effect similar to ones previously reported in biology (Bittihn et al., 2020;Sánchez-Taltavull et al., 2016).
Thus, we have determined the impact of a pandemic in the number of available healthcare workers for the different possible infection rates. It is important to note the different regimes. Low a and low b would correspond to an infectious disease that it is slowly spread among the healthcare workers, and few or no measures would be necessary. High a and high b would correspond to an infectious disease that spread rapidly and the desynchronization could help to keep the healthcare workforce available. High a and low b represent an epidemic with high case numbers among the general population but to which healthcare workers remain relatively well protected and avoid transmission amongst themselves. This would be the case in which high social distancing measures are applied in the hospital. Low a and high b would represent the opposite scenario -an epidemic that does not significantly reach the general population but spreads mainly in hospital settings -which has been seen in other outburst like MERS (Hastings et al., 2016).
An important limitation of our model, is that it assumes the infections coming from the general population to be constant. A possible future direction would be to build a more complex model in which the hospital dynamics are coupled to an epidemiological model of the general population (e.g. a city where the hospital is), and the pandemic in the general population and the hospital coevolves.
Recent models have revealed that the efficiency of desynchronization strategies is dependent on the frequency of rotations (Ely et al., 2020;Lim et al., 2020). Our study extends these models by incorporation assumptions on work efficiency and reinfection rates.
Our model has limitations that can be addressed in multiple ways. Another extension would be to add asymptomatic patients to the model, that is, patients that pass the infection without presenting symptoms. Furthermore, the model does neither address at what costs such a desynchronization strategy would come. Given that a complete set of workers is always present, there might be losses of acquired experience for the group absent. Furthermore, these models neither address the impact of length of duty shifts nor the impact of potentially impaired communication as a consequence of decreased interactions between healthcare workers.
In summary, our model is a starting point to study how to protect healthcare workers while determining economic impact during a pandemic outbreak.