Dynamic optimal allocation of medical resources: a case study of face masks during the first COVID-19 epidemic wave in the United States

1.   Introduction

At the beginning of a pandemic, the policymakers have to decide how to allocate limited medical resources among the general public and some special groups such as healthcare workers and vulnerable people. Together, mathematical models alongside a dynamic optimal control analysis will provide insights into finding the optimal allocation of medical resources. In this paper, we will use the COVID-19 pandemic ^[1] as an illustrative example to show how to solve the dynamic optimal allocation of face masks among healthcare workers and the general public.

During the first two months of 2020, the pandemic was mainly restricted to China ^[2,3]. The United States did not take very serious or effective actions until March 2020 ^[4], even though the first case was reported on January 22, 2020 ^[5]. Due to different cultures and medical systems, the control strategies varied from one country to another. For instance, wearing face masks in public areas is ubiquitous and even enforced by compulsory policies in China, Japan, and South Korea for everyone ^[6]. On the other hand, in March 2020, the US Surgeon General advised the general public not to buy face masks ^[7], and the Centers for Disease Control and Prevention (CDC) commented that face masks won't protect healthy people from getting the SARS-CoV-2 ^[8]. Nevertheless, according to the studies in ^[9,10], the aerosol transmission of SARS-CoV-2 is possible because the virus can last in aerosols for hours and remain viable and infectious. Moreover, there is plenty of scientific evidence showing that even though a person does not develop any symptom, they may still be infectious ^[11]. Thus, wearing face masks in public is an effective way to reduce the spreading of COVID-19, and more importantly, to prevent asymptomatic carriers ^[12] from infecting others.

Some dedicated compartmental models were developed in ^[13,14,15] to assess the community-wide impact of massively using face masks by the general public, which shows a high efficacy in curtailing community transmission and reduces the burden of the pandemic. In this paper, we will investigate how to optimize the allocation of face masks between healthcare workers (HCW) and the general public during the early stage of the COVID-19 outbreak. In April 2020, the CDC started to recommend wearing cloth face coverings for the general public (GP). Due to the initial supply shortage, the CDC advised saving face masks for healthcare workers ^[16]. It poses a natural question to ask whether one should distribute either a portion or none of the face masks to the general public during the epidemic outbreak. In other words, should the general public continue to reserve all the face masks solely for the HCW? We are particularly interested in the group of HCW since higher nurse-to-patient staffing ratios can result in healthier patients; hence, protecting the health of the limited number of HCW is paramount ^[17,18]. In general, it may take several years to train a qualified HCW to do a professional job. The objective of this work is to design an epidemic model with face masks and stay-at-home factors to find the optimal dynamical distribution of face masks among healthcare workers and the general public and to minimize the total number of deaths.

This paper is organized as follows. In Section 2, an optimal control model based on a new two-group SIR model is introduced, where the effects of wearing face masks and the stay-at-home policy are considered. The justification and estimation of the relevant parameters of our proposed model are discussed in Section 3. Section 4 presents the numerical results of the dynamic optimal control, where three different scenarios are demonstrated. Finally, some conclusions and discussion are given in Section 5.

2.   A two-group SIR model

Our model is based on the standard SIR (susceptible-infected-recovered) model proposed in ^[19]. To distinguish between two groups (i.e., the general public (GP) and healthcare workers (HCW)), we introduce six compartments $S_1$, $I_1$, $R_1$, $S_2$, $I_2$, and $R_2$ for the SIR classifications in the two groups. The subscripts $1$ and $2$ refer to the GP and the HCW, respectively. We do not include the asymptomatic (or latent) compartment because we assume that all infected individuals are infectious, even though they do not develop any symptoms ^[11]. Let $K_1(t)\ge 0$ and $K_2(t)\ge 0$ be the numbers of face masks distributed in the two groups over the given period $[0, T]$. The average number of available face masks for these two groups are

Here, we assume that the recovered GP ($R_1$) does not need face masks, but the recovered HCW ($R_2$) will return to their job and still need to wear face masks. Moreover, the infected HCW ($I_2$) are released from their duty of healthcare and thus considered as a part of the GP. Now, we can set the transmission rates for the GP and the HCW to be ${\beta} f_1(\rho_1)$ and ${\beta} f_2(\rho_2)$, where ${\beta} > 0$ is the intrinsic transmission rate of SARS-CoV-2, and $f_1$ and $f_2$ are positive and decreasing functions. For simplicity, we choose

where $1-{\alpha}\in(0, 1)$ stands for the maximum efficacy of wearing face masks, and $r > 1$ indicates the higher risk of infection for the HCW compared to the GP. Next, we introduce the ratio of infected individuals over the number of HCW:

and assume that the recovery and death rates depend on $\rho$ according to the two given functions ${\gamma}(\rho)$ and ${\delta}(\rho)$, respectively. It is reasonable to assume that the higher $\rho$ is, the smaller ${\gamma}(\rho)$ and the larger ${\delta}(\rho)$ will be. In our simulations, we will choose

where ${\gamma}_0 > 0$ and ${\delta}_0 > 0$ are the recovery and death rates of an infected individual with sufficient healthcare ($\rho\to0$), respectively, and ${\gamma}_\infty > 0$ and ${\delta}_\infty$ are the rates when the healthcare system is overwhelmed ($\rho\to\infty$). Obviously, ${\gamma}_0 > {\gamma}_\infty$ and ${\delta}_0 < {\delta}_\infty$. To estimate the efficacy of the stay-at-home policy, we introduce another parameter $q > 0$ to account for the portion of susceptible individuals who limit their activity in an isolated region (staying at home for example) away from the infected group. Finally, we are ready to state our two-group SIR model which consists of six ordinary differential equations

where $\rho_1$, $\rho_2$, and $\rho$ are defined in (2.1) and (2.3), respectively, the functions $f_1$ and $f_2$ are given in (2.2), and the functions ${\gamma}$ and ${\delta}$ are defined in (2.4). The prime symbol on the left-hand side of the equations denotes the ordinary derivative in time $t$. Such two-group SIR epidemic models were widely known in the literature, see e.g., ^[20] and references therein. Note that $f_i(\rho)$ (with $i = 1$ or $i = 2$) appears twice in the incidence rates. This is because face masks are more effective in reducing disease transmission if both susceptible and infected individuals are wearing them.

Let $T$ be the considered duration of an outbreak and $K_{\max}$ be the maximum capacity of the daily production number of face masks. The objective is to minimize the total number of deaths

subject to the point-wise control constraints describing the production capacity limits

The above optimization model gives a nonlinear constrained optimal control problem, whose optimal control (may not be unique) can be mathematically characterized via Pontryagin's minimum principle ^[21]. Due to its high non-linearity, we will numerically solve the above optimal control problem using a direct transcription method ^[22,23] based on Legendre-Gauss-Radau pseudo-spectral collocation ^[24,25] and nonlinear programming (NLP) solvers (e.g., IPOPT ^[26] and SNOPT ^[27]). Different from indirect methods that require to derive necessary optimality conditions, such direct transcription methods are more flexible in treating extra constraints and are widely supported by well-developed general optimal control software packages that are ready to be used. In particular, our following numerical simulations were performed with the free and open-source Imperial College London Optimal Control Software (ICLOCS2) ^[28] within MATLAB, where the derivatives were numerically computed by the algorithmic differentiation toolbox Adigator ^[29].

3.   Parameter estimation

There have been plenty of works dedicated to the estimation of key epidemic parameters such as the basic reproduction number and the serial interval during the early outbreak ^[30,31,32]. To estimate the values of the parameters in our proposed model, we shall fit the reported data on cumulative confirmed case numbers, denoted by $C(t)$. The change rate of $C(t)$ is the same as the newly infected case number, which according to our model is

The initial numbers of the susceptible population and healthcare workers are estimated as $S_1(0) = 310, 000, 000$ ^[33] and $S_2(0) = 16, 000, 000$ ^[34]. Throughout this paper, the time unit is in days. We chose March 13, 2020 as the initial time $t = 0$ since an emergency declaration was warranted for the COVID-19 pandemic on that date ^[4]. Accordingly, we set $I_1(0) = 1896$ ^[5]. Since the epidemic wave started and stay-at-home policy was adopted following the declaration, we let $I_2(0) = R_1(0) = R_2(0) = 0$. Based on the cumulative case numbers reported from February 28, 2020 to March 13, 2020 ^[5], we can estimate the intrinsic growth rate as ${\beta} S_1(0)-{\gamma}_0-{\delta}_0 = 0.35$ by fitting a simple exponential growth model

which can be obtained by approximating our model under the assumptions that, before March 13, 2020, the stay-at-home policy was not adopted ($q = 0$), few people wore face masks ($K_1+K_2\ll S_2\ll S_1$), and the infected population was small ($I_2\ll I_1\ll S_2$). Consequently, we have $S_1(t)\approx S_1(0)$ and $\rho_1, \rho_2, \rho, I_2\approx0$. The equation (2.7) is approximated by a linear equation $I_1' = [{\beta} S_1(0)-{\gamma}_0-{\delta}_0]I_1$, with the solution $I_1(t) = I_1(0)e^{[{\beta} S_1(0)-{\gamma}_0-{\delta}_0]t}$. The equation for the cumulative case numbers $C(t)$ in (3.1) can be approximated as $C'(t) = {\beta} S_1(0)I_1(t) = {\beta} S_1(0)I_1(0)e^{[{\beta} S_1(0)-{\gamma}_0-{\delta}_0]t}$. This, together with the initial condition $C(-\infty) = 0$, gives (3.2).

The basic reproduction number $R_0 = {\beta} S_1(0)/({\gamma}_0+{\delta}_0)$ for SARS-CoV-2 was estimated to be approximately $2.2$ ^[35]. In view of ${\beta} S_1(0)-{\gamma}_0-{\delta}_0 = 0.35$, we obtain ${\gamma}_0+{\delta}_0 = 0.29$ and ${\beta} S_1(0) = 0.64$ from a simple calculation. To estimate ${\gamma}_0$ and ${\delta}_0$, we need to use the death rate, which varies from $4%$ to $7.5%$ ^[36]. Here, we chose ${\delta}_0/({\gamma}_0+{\delta}_0) = 7%$. Consequently, we have ${\delta}_0 = 0.02$ and ${\gamma}_0 = 0.27$. There is no data available for the case when the healthcare system is overwhelmed by too many infected patients, and we simply assume that with an overwhelmed healthcare system, the death rate will reach ${\delta}_\infty = 0.1$ and the recovery rate reduces to ${\gamma}_\infty = 0.1$.

The values for ${\alpha}$ and $r$ cannot be found in the literature. Here, we set ${\alpha} = 0.9$ (face masks can at most reduce the transmission rate by $10%$) and $r = 3$ (the healthcare workers are three times more likely to be infected than the general public). We have used some other values for these parameters, and the results do not vary much. Especially, we numerically find that the optimal distributions of face masks have very similar patterns for other reasonable choices of ${\alpha}$ and $r$.

To estimate the efficacy of stay-at-home parameter $q$, we assume $\rho\approx0$ and ignore the population of healthcare workers ($S_2 < < S_1$ and $I_2 < < I_1$). For simplicity, we also assume that the majority do not wear face masks ($\rho_1\approx0$) and then fit the confirmed cumulative case numbers reported from March 13, 2020 to April 30, 2020 ^[5] by the following simplified SIR model:

where $S(0) = 310, 000, 000$, $C(0) = 1896$, ${\beta} = 0.64/S(0)$, ${\gamma}_0 = 0.27$, and ${\delta}_0 = 0.02$. The equation for $R$ can be decoupled from the system. The two unknown parameters $q$ and $I(0)$ are estimated via MATLAB's nonlinear least-squares curve fitting solver $\texttt{lsqnonlin}$. With a random initial guess, the $\texttt{lsqnonlin}$ solver converges to the (local) optimally estimated parameters: $q = 0.03$ and $I(0) = 1070$. Numerically, we observe the fitted parameters $q$ and $I(0)$ are insensitive to the chosen initial guess, as well as the different possible choices of $R(0)$. It is tempted to treat $R(0)$ as an additional unknown parameter, but in this way, its fitted value is not uniquely determined by the given data. Hence, we reasonably set $R(0) = 0$ by assuming nobody recovered at the beginning. Figure 1 illustrates the close match between the reported $C(t)$ data and the simulated $C(t)$ based on the fitted parameters. We point out that the stay-at-home term $qS$ is crucial for achieving such a satisfactory fitting while fixing the other model parameters. The assumed and estimated parameter values are listed in Table 1.

4.   Dynamic optimal control

In this section, we report some inspiring simulation results based on our proposed optimal control model for optimizing the allocation of face masks among HCW and the GP with our estimated parameters. We highlight that in the application of our model to different countries/areas, one may need to re-estimate some of the parameters based on the reported data and the population scale. Our tested choices of the maximum daily production capacity of the face masks are only for demonstrating our proposed model, which does not reflect real-life situations. Moreover, our current model does not take those homemade cloth face masks and the effects of mandatory quarantine into account. Therefore, our model outcomes are mainly for a qualitative comparative analysis.

We set $T = 100$ and consider 3 scenarios of the maximum daily production capacity of face masks:

(i) $K_{\max} = 16, 000, 000$ (each HCW can have at most one face mask every day);

(ii) $K_{\max} = 80, 000, 000$ (each HCW can have at most five face masks every day);

(iii) $K_{\max} = 160, 000, 000$ (each HCW can have at most ten face masks every day).

For scenario (i) with a very limited supply, our optimal control results are reported in Figure 2. It shows that the face masks should be distributed to the HCW only, which agrees with the guideline from the CDC. With this policy control, the total death numbers of the GP and the HCW are $J_1 = 386, 017$ and $J_2 = 104, 707$, respectively. With more investments from the government and industries, the supply of face masks have been quickly boosted to a higher level. In this context, we would want to ask whether this CDC guideline continues to be optimal with an increased supply of face masks. The simple answer is no, as clearly shown in the next two scenarios.

For scenario (ii) with moderate supply, our optimal control results are reported in Figure 3. Quite different from scenario (i), it shows that the face masks should be approximately equally distributed at the beginning, and then gradually shifted to the HCW, though not all face masks are given to the HCW across the whole pandemic outbreak. The total death numbers of the GP and the HCW are $J_1 = 144, 818$ and $J_2 = 27, 145$, respectively. As a comparison, Figure 4 shows the corresponding outcomes if strictly following the CDC guidelines to allocate all face masks to HCW, where the total death numbers of the GP and the HCW become $J_1 = 182, 409$ (with a 26% increase) and $J_2 = 29, 890$ (with a 10% increase), respectively. The significant increase (about 23%) in the total death toll is somewhat surprising but reasonable, since an increased number of infected GP due to not wearing face masks would further infect more HCW, and ultimately lead to higher death rates and more deaths in both groups.

For scenario (iii) with very sufficient supply, our optimal control results are reported in Figure 5. Different from both scenarios (i) and (ii), it shows that the majority of face masks should be distributed to the GP at the beginning of the epidemic outbreak, and then gradually shifted to the HCW during the outbreak. The total death numbers of the GP and the HCW are $J_1 = 80, 897$ and $J_2 = 12, 208$, respectively. For comparison, Figure 6 shows the corresponding outcomes if following the CDC guideline to allocate all face masks to HCW, where the total death numbers of the GP and the HCW become $J_1 = 155, 540$ (with a 92% increase) and $J_2 = 21, 871$ (with a 80% increase), respectively. Astonishingly, the total death numbers almost doubled if strictly following the reasonably sound CDC guidelines, where, based on our model setting, the protective effect of wearing too many face masks for HCW is essentially saturated. It is also worthwhile to notice that the optimal allocation of masks has greatly flatted the infectious peak curves and hence lead to much fewer deaths.

As an alternative illustration, the relevant functions $f_1$, $f_2$, $\gamma$, $\delta$ with respect to time are compared in Figure 7, where we indeed observe slightly higher death rates and lower recovering rates if all face masks are reserved solely for HCW (as advised by CDC). Scenario (iii) indicates that an appropriately balanced allocation of face masks between the GP and the HCW plays a significant role in saving more lives. In summary, the three different scenarios manifest that the allocation of face masks needs to be carefully optimized, especially when the supply becomes gradually more sufficient.

5.   Conclusions and discussion

In this paper, we constructed a two-group SIR model to optimize the distribution of face masks among healthcare workers (HCW) and the general public (GP). When the supply of face masks is short, our results indicate that all face masks should be reserved for the HCW. This coincides with the advice from the CDC in April 2020 ^[16]. However, when there is a sufficient supply of face masks, the general public should share a large portion of face masks at the beginning of an epidemic outbreak. This result somewhat contradicts the recommendations given in March 2020 by the US Surgeon General ^[7] and the CDC ^[8]. The optimality of this reasonable and sound CDC guideline highly depends on the supply level of face masks, which changes frequently and varies by location, and hence this guideline should be modified according to the supply of face masks. Based on our choices and estimations of parameter values, assuming that the supply of face masks is sufficient, and the stay-at-home policy remains effective, our model indicates that the first epidemic wave would have ended in May 2020, with a cumulative total 93, 105 deaths (80, 897 deaths from the GP and 13, 208 deaths from the HCW). Note that the stay-at-home policy was released before the end of the first epidemic wave, and an even stronger second epidemic wave arose afterward. Unlike physical phenomena which can be observed from repeated experiments, epidemic outbreaks cannot be tested multiple times. Thus, it is important to explore the effects of hypothetical control measures on epidemic waves ^[37]. Our model analysis and numerical simulation provides theoretical experiments on what would have occurred if the general public were advised to wear face masks at the beginning of the first epidemic wave and the stay-at-home policy was enforced until the end of the epidemic wave.

There are some limitations to our studies. For instance, the values for parameters ${\delta}_\infty$, ${\gamma}_\infty$, ${\alpha}$, and $r$ are arbitrarily chosen. More real-life data collections are required to obtain more reliable estimations on the death and recovery rates of SARS-CoV-2 in the case of the overwhelmed healthcare system, the efficacy of wearing face masks, and the risk to healthcare workers. With these limitations being said, our model predictions may vary a lot whenever the involved parameters are changed. In particular, the reopening of the country will significantly diminish the stay-at-home efficacy and generate multiple epidemic waves.

Our proposed model can be extended in several different ways. For example, it is interesting to alternatively enforce an integral constraint $\int_{0}^t K_1(\tau)+K_2(\tau) d\tau \le (K_{0}+K_{\max}t),$ which is more practical since it allows for the flexible usage of face masks according to the epidemic dynamics and possible initially stocked face masks (denoted by $K_0$). Our developed model can be generalized to optimize the allocations of other personal protective equipment (PPE) among various groups of populations (e.g., public workers vs the general public). It is also possible to further differentiate the protection effects between different types of face masks. In addition, the objective function may also include some economic considerations, such as the production and transportation costs of face masks. The generalization of our model to optimally allocate the limited vaccination is also interesting.

Acknowledgments

We are very grateful to the two referees for their careful reading and valuable suggestions which have helped to improve the presentation of this paper. X. S. Wang acknowledges the support from the Louisiana Board of Regents Support Fund under contract No. LEQSF(2022-25)-RD-A-26.

Conflict of interest

Xiang-Sheng Wang is a guest editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. All authors declare no conflicts of interest in this paper.