Balancing mitigation strategies for viral outbreaks

1.   Introduction

Control and elimination of infectious diseases have been a focus of public health officials since the early 1950s. As antibiotics ^[1,2], sanitation ^[3], and vaccinations ^[4,5,6] were introduced, eradicating diseases became more feasible ^[7,8]. Nevertheless, several factors have led to the emergence of new infectious diseases and the reemergence of existing ones, including resistance of microorganisms to medication ^[9,10,11], demographic evolution ^[12], urbanization ^[13,14], and increased travel ^[15]. These diseases include Lyme disease in 1975 ^[16], Legionnaires disease in 1976 ^[17], toxic shock syndrome in 1978 ^[18], Hepatitis C in 1989 ^[19], Hepatitis E in 1990 ^[20], and Hantavirus in 1993 ^[21]. HIV (human immunodeficiency virus) emerged in 1981 as one of the most alarming sexually transmitted diseases in the world ^[22,23]. As the result of antibiotic resistance, tuberculosis, pneumonia, and gonorrhea are reemerging ^[24]. Because of climate change, malaria, dengue, and yellow fever have also reemerged and are spreading into new areas ^[25]. Moreover, it is not unusual for diseases such as plague and cholera to erupt from time to time ^[26]. The reemergence of the Ebola virus disease (EVD) in 2013 has puzzled the world ^[27]. Fifteen million deaths are directly related to the reemergence of infectious diseases every year, which remains a serious medical burden around the globe ^[28,29].

In the field of epidemiology, the application of optimal control strategies has proven to be a critical tool in managing and mitigating the spread of infectious diseases ^[30,31,32]. Optimal control theory provides a robust framework that allows one to incorporate one or more time-dependent control functions into a nonlinear dynamic system to achieve the best possible outcome for a specified objective ^[33,34]. It uses mathematical and computational techniques aimed at identifying the most effective interventions, such as vaccination ^[35], treatment plans ^[36], quarantines, and social distancing ^[37], while optimizing the use of resources ^[38,39]. By integrating optimal control into disease models, researchers can predict and influence the progression of an outbreak, ultimately reducing morbidity, mortality, and economic impact ^[40,41,42].

Optimal control problems have been extensively studied across various infectious diseases. For instance, during the 2014–2016 Ebola outbreak in West Africa, optimal control models were developed to evaluate the effectiveness of different intervention strategies, including isolation ^[43], contact tracing ^[44], and public health campaigns ^[45,46]. Dengue fever, a mosquito-borne viral disease, also presents a significant public health challenge in tropical and subtropical regions ^[47]. Optimal control strategies for dengue have focused on reducing mosquito populations and limiting human exposure through targeted insecticide use ^[48], environmental management ^[49], and vaccination ^[50,51]. The COVID-19 pandemic has brought the application of optimal control strategies to the forefront of public health efforts globally ^[34,52]. Researchers have developed models to determine the best combination of non-pharmaceutical interventions (NPIs) such as social distancing ^[53], lockdowns ^[54], and mask-wearing ^[55], alongside vaccination rollouts ^[54,56]. These models have been instrumental in guiding policy decisions, helping to balance the trade-offs between controlling the virus and minimizing societal disruption ^[57,58,59]. Other diseases, such as influenza ^[60,61], malaria ^[62], Zika virus ^[63], various types of cancer, and HIV ^[64], have also been the focus of optimal control studies, each contributing to the development of a robust framework for disease management ^[65,66].

In this paper, we develop a robust optimal control algorithm to simulate five epidemic scenarios: "no control", "reconstructed", "social distancing control", "vaccine control", and "both controls concurrently", using biological models for the post-vaccination stage of a viral outbreak. For numerical validation, we employ real data for the SARS-CoV-2 Delta variant in different regions of the United States of America from July 9, 2021 to November 25, 2021. By comparing several control scenarios, we provide an overarching study of crucial intervention strategies, offering valuable insights into disease transmission. Our innovative approach to modeling the cost of control gives rise to a fast trust-region optimization procedure for solving a broad range of nonlinear control problems.

The paper is organized as follows: Section 2 outlines mathematical preliminaries essential to our analysis. We introduce and examine the social distancing control strategy in Section 3. In Section 4, the vaccine control is investigated, and Section 5 is dedicated to the combined implementation of social distancing and vaccination controls. In Section 6, we compare the aforementioned scenarios and assess their efficiency. Additional figures and tables are presented in the Appendix.

2.   Mathematical preliminaries

Mathematical models have been extensively used to analyze the spread of infectious diseases, such as Ebola, dengue, plague, influenza, and COVID-19, and to forecast their impact ^[67,68]. Disease models range in complexity from basic Susceptible - Exposed - Infectious - Recovered (SEIR) ODE systems with few parameters to more sophisticated models that include numerous compartments and parameters to account for isolation ^[69], hospitalization, testing, contact tracing ^[70], social distancing ^[71,72], vaccination ^[73], and others ^[74,75]. In this study, we consider the $\hat{S}_u\hat{S}_v\hat{I}_u\hat{I}_v\hat{R}\hat{D}$ compartmental model that was first introduced in ^[76] and further studied in ^[77]:

This model incorporates the vaccination status of both susceptible and infected individuals, as well as the possibility of losing immunity and becoming infected among both vaccinated and unvaccinated populations ^[76]. Epidemic model (2.1) has 6 compartments: Susceptible unvaccinated ($\hat{S}_u$), Susceptible vaccinated ($\hat{S}_v$), Infected unvaccinated ($\hat{I}_u$), Infected vaccinated ($\hat{I}_v$), Recovered ($\hat{R}$), and Deceased ($\hat{D}$). In (2.1), the time-dependent parameter, $\zeta(t)$, is the disease transmission rate, $\mu$ is the average vaccination rate, $\delta _r$ is the rate at which individuals lose immunity after acquiring the virus, $\delta_v$ is the waning vaccine immunity rate, and $\gamma_{u, r}$ and $\gamma_{v, r}$ are the recovery rates for unvaccinated and vaccinated individuals, respectively. The virus death rates for unvaccinated and vaccinated humans are denoted by $\gamma_{u, d}$ and $\gamma_{v, d}$, and $\alpha$ ($0 < \alpha < 1$) is the measure of vaccine efficacy for the current strain. In the above, $t\in[0, T]$, where $T$ is the length of the study period.

In model (2.1), it is assumed that the duration of each virus strain is relatively short compared to the time it takes for the population to change due to birth, migration, death of causes rather than the virus, etc. Therefore natural birth and death are omitted in the model, and we assume that at any given time, $t$, the population of the region is $N-D(t)$, where $N$ is the population at $t = 0$, that is, $N = \hat{S}_u(0)+\hat{S}_v(0)+\hat{I}_u(0)+\hat{I}_v(0)+\hat{R}(0)+\hat{D}(0)$. We also assume that vaccination is applied only to susceptible individuals, and infected and recovered individuals are not vaccinated until they lose immunity and move back to the susceptible class. So, this model does not account for asymptomatic infected individuals that can be vaccinated while still infected or in the recovered stage. According to model (2.1), the loss of immunity after infection is the same for both vaccinated and unvaccinated individuals.

In ODE system (2.1), $\zeta(t)$ is a time-dependent parameter that accounts for real-life nonmedical preventive measures (social distancing, bans on travel and large gatherings, handwashing, etc.) aimed at bringing down the disease transmission. In our numerical simulations, $\zeta(t)$ is reconstructed for individual states from real data on new incidence cases and daily deaths ^[78] using the algorithm proposed in ^[77]. Another important parameter, $\mu$, is the average vaccination rate, which is pre-estimated based on CDC reports ^[76,79] (for individual states) by dividing the change in the percentage of vaccinated people at the start and at the end of the study window by the length of this window. As proposed in ^[76,80,81], we assume that susceptible unvaccinated individuals move to the susceptible vaccinated class at the rate proportional to the current number of susceptible unvaccinated individuals. By estimating $\zeta(t)$ and $\mu$ in this manner, we analyze the actual disease progression (what is called the "reconstructed" scenario in our experiments).

The goal of the optimal control problem is to see how much the real-life scenario can be improved through the implementation of control strategies that optimize a specific objective functional. To that end, in the "real-life" model (2.1), we replace $\zeta(t)$ with $\beta (1-u_1(t))$, where $\beta$ is the original disease transmission rate and $u_1(t)$ is the social distancing control that we aspire to optimize. Clearly, one has $\beta\ge \max_{t\in [0, T]}\zeta(t)$ since $\zeta(t) = \beta(1-\tilde u_1(t))$, where $\tilde u_1(t)$ is a "real-life" non-optimal social distancing control ($0\le \tilde u_1(t) < 1$ for $t\in [0, T]$). In order to facilitate the vaccination strategy, we replace $\mu$ with $\nu u_2(t)$, where $\nu > \mu$ is the capacity of vaccination and $u_2(t)$ is the vaccination control. The feasible set for $u_1(t)$ and $u_2(t)$ is

To simplify the biological model, we normalize the state variables, $S_u : = \frac{\hat{S}_u}{N},$ $S_v : = \frac{\hat{S}_v}{N},$ $I_u : = \frac{\hat{I}_u}{N},$ $I_v : = \frac{\hat{I}_v}{N},$ $R: = \frac{\hat{R}}{N},$ and $D : = \frac{\hat{D}}{N},$ and arrive at the following controlled system of differential equations:

The two control functions, $u_1(t)$ and $u_2(t)$, are intended to lower the normalized force of infection, $\beta(1-u_1(t))\frac{S_u(t)+(1-\alpha)S_v(t)}{1-D(t)}(I_u(t)+I_v(t))$, while keeping the costs at bay. The costs of control are understood in a general sense, which includes a negative impact on the economy, mental health, education, and other aspects of life ^[82]. With that in mind, we propose the following objective functional

where $\mathbf{x}(t): = [S_u(t), S_v(t), I_u(t), I_v(t), R(t), D(t)]^\top$ combines all normalized state variables of the model, $\boldsymbol u(t) : = [u_1(t), u_2(t)]^\top$ represents the two control strategies, $\mathbf{c(u)} : = [c_1(u_1), c_2(u_2)]^\top$ stands for the vector of cost functions associated with $u_1(t)$ and $u_2(t)$, respectively, and $\boldsymbol \lambda : = [\lambda_1, \lambda_2]^\top$ is the weight vector for the cost functions, $c_i(u_i)$, $i = 1, 2$. Thus, the optimal control problem for what we call the "both controls concurrently" scenario is to minimize objective functional (2.4) subject to ODE system (2.3).

In the next two sections, we will also introduce two special cases, "social distancing control" (i.e., social distancing control only) and "vaccine control" (i.e., vaccine control only), with functional (2.4) and model (2.3) adjusted accordingly. Numerical results for these three control problems will be compared with the aforementioned "reconstructed" scenario (2.1) and with the "no control" scenario, described by $S_u I_u R D$ model (2.5) below. To introduce the "no control" disease progression, we assume that neither social distancing nor vaccination control is enforced, that is, the epidemic is running its course. Hence, the disease transmission rate is constant (and equal to $\beta$) and there are no vaccinated compartments for either susceptible or infected humans:

According to Pontryagin's Minimum Principle ^[33,83], if ${\mathbf u} \in \mathcal{U}$ is an optimal control strategy with respect to the objective functional $J({\mathbf x}, {\mathbf u}) = h({\mathbf x}(T)) + \int_{0}^{T} L({\mathbf x}(t), {\mathbf u}(t))\, dt$ and the system of equations $\dot{\mathbf x} = f({\mathbf x}, {\mathbf u})$, ${\mathbf x}(0) = {\mathbf x_0}$, then there is a trajectory ${\mathbf p}(t)$ such that

Therefore, solving the optimal control problems for the "social distancing control", "vaccine control", and "both controls concurrently" scenarios comes down to minimizing the Hamiltonian, $H({\mathbf x}, {\mathbf u}, {\mathbf p})$, with respect to ${\mathbf u}$ subject to costate system (2.6) and the biological model $\dot{\mathbf x} = f({\mathbf x}, {\mathbf u})$, ${\mathbf x}(0) = {\mathbf x_0}$. The complexity of this minimization problem largely depends on the choice of the cost function, $\mathbf{c(u)}$, in the objective functional (2.4). In our algorithms, we employ a twice continuously differentiable cost function, $\mathbf{c(u)} : = [c_1(u_1), c_2(u_2)]^\top$, with the following key properties ^[82]:

● $c_i(u)$, $i = 1, 2$, are defined in $\mathcal{D}_i\supseteq [0, 1)$ ensuring that each domain contains the feasible set;

● $c_i(0) = 0$, $i = 1, 2$, implying zero cost when no control is applied;

● $\lim_{u \rightarrow 1^{-}} c_i(u) = \infty$, guaranteeing that the cost becomes prohibitive as $u$ approaches the upper bound of the control range;

● $c_i^{\prime}(u) > 0$ when $u > 0$, and $c_i^{\prime}(u) < 0$ when $u < 0$, suggesting that $c_i(u)$ is increasing for positive control values and preventing $u$ from becoming negative;

● $c_i^{\prime\prime}(u) > 0$ for all $u \in \mathcal{D}_i$, indicating that $c_i(u)$ is strictly convex, $i = 1, 2$.

Note that $c_i(u)$, $i = 1, 2$, satisfying the above conditions, are necessarily nonnegative in their respective domains, $\mathcal{D}_i$. In the numerical simulations presented in this paper, for different control scenarios, we consider the following four cost functions, $i = 1, 2$ ^[84]:

It is important to mention that for $c_{i, 2}(u) = u^2$ ^[85,86], which is often used in optimal control problems (or for a more general function, $c_{i, 2}(u) = w_1u+w_2u^2$ ^[87,88]), the requirement $\lim_{u \rightarrow 1^{-}} c_i(u) = \infty$ is not met. In our experimental framework, $c_{i, 2}(u) = u^2$, $i = 1, 2$, is used for comparison to highlight the importance of the property $\lim_{u \rightarrow 1^{-}} c_i(u) = \infty$. For $c_{i, 1}$, $c_{i, 3}$, and $c_{i, 4},$ all assumptions listed above are fulfilled and, as our experiments show, all candidates for the global minimum of $H({\mathbf x}, {\mathbf u}, {\mathbf p})$ with respect to ${\mathbf u}$ subject to the state and costate systems are feasible, i.e., inequality constraints, $u_i\ge 0$ and $u_i < 1$ (or control-specific constraint, $u_i < b_i$ ^[87,88]), do not have to be enforced in the optimization algorithm. This is not the case with $c_{i, 2}(u) = u^2$. For this cost function, one or both coordinates of the global minimum are often greater than $1$, especially for small values of $\lambda_i$. The details are presented in Sections 3–5.

For the three optimal control problems "social distancing control", "vaccine control", and "both controls concurrently", investigated in Sections 3–5, respectively, we use $c_{1, j}(u) = c_{2, j}(u)$, $j = 1, 2, 3, 4$. However, in some applications, the cost functions associated with different controls may need to be different. The computational algorithm can easily be adapted for that. The graphs of the cost functions, $c_{i, j}(u)$, $u\in [-1, 1)$, $i = 1, 2$, $j = 1, 2, 3, 4$, are shown in Figure 1.

3.   Social distancing control

In this section, we consider an intervention scenario where only social distancing controls are implemented. Similar to (2.5), there are no vaccinated compartments for either susceptible or infected humans, and the model with social distancing control, $u_1(t)\in \mathcal{U}$, takes the form:

A higher control value implies more strict social distancing measures, leading to a reduced transmission. On the other hand, $u_1(t) = 0$ yields no social distancing and a disease transmission rate reaching its full potential. With no vaccination control, the normalized force of infection is $\beta(1-u_1(t))\frac{S_u(t)I_u(t)}{1-D(t)}$, and one arrives at the following objective functional:

designed to achieve a balance between reducing the spread of the disease and managing the associated mitigation costs over the period \([0, T]\). In (3.2), $\mathbf{x}: = [S_u, I_u, R, D]^\top$. Equations (2.7), (3.1), and (3.2) give rise to the Hamiltonian:

where $\mathbf{p}: = [p_1, p_2, p_3, p_4]^\top$. By Pontryagin's Minimum Principle ^[33,83], $u_1 = \operatorname*{arg\, min}_{v \in \mathcal{U}} H({\mathbf x}, v, {\mathbf p}),$ subject to state system (3.1) and costate system (2.6) in the form

This leads to the following $2^{\text{nd}}$-order numerical algorithm for nonlinear constrained minimization:

In Algorithm 1, $F(\theta)$ is a discrete analog of the partial derivative of the Hamiltonian, $H({\mathbf x}, u_1, {\mathbf p})$, with respect to $u_1$, $J(\theta)$ is the Jacobian of $F(\theta)$, $I$ is the identity matrix in the solution space, $\varrho$ is the step size, and $\omega_k$ is the regularization sequence. The derivative of $H({\mathbf x}, u_1, {\mathbf p})$ with respect to $u_1$ exists, since $c_{1, j}(u)$, $j = 1, 2, 3, 4$, are twice continuously differentiable by our assumption. In all our experiments, shifted Legendre polynomials were used to project the control function, $u_1(t)$, onto a finite dimensional subspace with $\theta$ being a vector of expansion coefficients. MATLAB's built-in function "ode15s" was employed to solve both ODE systems, (3.1) and (3.4), while "lsqnonlin" implemented the Levenberg-Marquardt optimization procedure.

To illustrate the efficiency of various control strategies, we use real data for the SARS-CoV-2 Delta variant of the COVID-19 pandemic in Alabama and Maryland from July 9, 2021 to November 25, 2021 ^[78]. In this section, the "social distancing control" scenario is compared to what we call "reconstructed" and "no control" scenarios. The "reconstructed" (or "real-life") scenario is described by system (2.1), where pre-estimated parameter values for the state of Alabama are set at $N = 5,031,362$, $\gamma_{u, r} = (1-0.005)/10$, $\gamma_{u, d} = 0.005/18.5$, $\gamma_{v, r} = (1-0.005/12.7)/10$, $\gamma_{v, d} = 0.005/18.5/12.7$, $\delta_r = 1/90$, $\delta_v = 0$, $\alpha = 0.8$, and $\mu = 0.0009143$ ^[76,77], while the time-dependent transmission rate, $\zeta(t)$, and case reporting rate, $\psi$, are reconstructed from CDC data ^[78] on daily new infections and deaths by regularized optimization algorithm. The reconstructed value of $\psi$ is equal to 0.154 (95%CI:[0.149, 0.159]) ^[77]. The initial values for (2.1) are $S_u(0) = 3,402,668/N$, $S_v(0) = 1,626,323/N$, $I_u(0) = 1584/N$, $I_v(0) = 787/N$, $R(0) = 0$, and $D(0) = 0$.

In the "reconstructed" or "real-life" scenario we assume that both social distancing and vaccination controls are present, but their implementation is not optimal and mimics real-life interventions put in place from 7/9/2021 to 11/25/2021. The hypothetical "no control" scenario is given by model (2.5), where neither social distancing nor vaccination control is applied. The hypothetical "social distancing control" scenario (3.1) represents the case where optimal social distancing control is implemented with no vaccination available (see Algorithm 1). The initial values for the state variables in these two cases are $S_u(0) = (3,402,668+1,626,323)/N$, $I_u(0) = (1584+787)/N$, $R(0) = 0$, and $D(0) = 0$. In (2.5) and (3.1), the constant transmission rate, $\beta$, is set to $0.416 = \max_{t\in [0, T]}\zeta(t)$ since $\zeta(t) = \beta(1-\tilde u_1(t))$, where $\tilde u_1(t)$ is a "real-life" non-optimal social distancing control ($0\le \tilde u_1(t) < 1$ for $t\in [0, T]$), as mentioned in Section 2 above.

Figure 2 shows that without any control measures in place, the daily number of infected people would be quite alarming until the strain runs its natural course. The hypothetical number of deceased individuals in the "no control" environment is dangerously high. On the other hand, with optimally enforced social distancing, even without vaccination, the daily number of infected and deceased humans is very low (though some cases appear to be delayed rather than prevented). The figure underscores the importance of social distancing in mitigating the impact of infectious disease outbreaks. The "reconstructed" curves in Figure 2 illustrate that the real-life control measures in the state of Alabama, which included both vaccination and social distancing, were very effective and saved many lives.

The first graph in Figure 3 demonstrates that optimal social distancing must be strictly enforced at the early ascending stage of a new wave in order to contain the virus. However, toward the end of the study period, the intensity of optimal control goes down due to its negative impact on the economy and overall quality of life. The figure shows that, for this particular value of $\lambda_1 = 0.01$, the control strategy for all four cost functions (2.8) remains within the feasible set the entire time, without inequality constraints enforced. At the same time, according to the second graph in Figure 3, the control strategy associated with $c_{1, 4}$ may not be a global minimum. One can see in Figure 1 that the cost of the $4^{th}$ control is the highest among all costs considered. Therefore, superior results for the $4^{th}$ control, shown in Figure 2, may be at the expense of the algorithm not sufficiently reducing the negative impact.

Figure 4 compares the actual number of new daily infections and deaths in the state of Alabama from July 9, 2021 to November 25, 2021 ^[78], to new infections and deaths in the case of hypothetical optimally controlled social distancing (for the four cost functions (2.8)). It also illustrates the model fit with 100 bootstrapping iterations for uncertainty quantification ^[76,77]. The figure supports our earlier observation that optimal implementation of social distancing prevents a considerable number of deaths. It illustrates that the daily number of newly infected people at the early stage of the cycle, where the optimal social distancing control is strictly enforced, is close to zero. However, new incidence cases and deaths increase toward the end of the study period, when the intensity of optimal control goes down. The surge of new infections and deaths in the last 40 days of the interval underlines the importance of a vaccination campaign for the prevention of cases. Optimal social distancing control is very powerful initially, but it is not sustainable for a long time.

Figures 5–7 represent the "no control", the "reconstructed", and the optimal "social distancing control" scenarios for the state of Maryland during the COVID-19 Delta variant from July 9, 2021 to November 25, 2021. By comparing Figures 2–4 to Figures 5–7, we can examine how regional differences, such as varying population density, healthcare infrastructure, and socioeconomic factors, influence the optimal control strategy ^[89,90].

For the state of Maryland, the "reconstructed" scenario, described by system (2.1), illustrates the efficiency of real-life interventions. In Maryland model (2.1), some parameters are the same as in the case of Alabama and others are different (i.e., state-specific). The pre-estimated parameter values for the state of Maryland are $N = 6,173,205$, $\delta_r = 1/90$, $\gamma_{u, r} = (1-0.005)/10$, $\gamma_{u, d} = 0.005/18.5$, $\gamma_{v, r} = (1-0.005/12.7)/10$, $\gamma_{v, d} = 0.005/18.5/12.7$, $\delta_v = 0$, and $\mu = 0.0007286$ ^[76,77]. The time-dependent transmission rate, $\zeta(t)$, and case reporting rate, $\psi$, are reconstructed from CDC data ^[78] on daily new infections and deaths by a regularized optimization algorithm. The reconstructed value of $\psi$ for Maryland is equal to 0.182 (95%CI:[0.172, 0.192]), see ^[77]. Pre-estimated initial values for the coordinates of $\mathbf{x}: = [S_u, S_v, I_u, I_v, R, D]^\top$ in (2.1) are $S_u(0) = 2,727,503/N$, $S_v(0) = 3,445,221/N$, $I_u(0) = 207/N$, $I_v(0) = 274/N$, $R(0) = 0$, and $D(0) = 0$.

As before, the hypothetical "no control" scenario is given by model (2.5) with no control applied, and the hypothetical "social distancing control" scenario (3.1) is the case where optimal social distancing control is employed with no vaccination available (see Algorithm 1). The initial values for systems (2.5) and (3.1) are $S_u(0) = (3,402,668+1,626,323)/N$, $I_u(0) = (1584+787)/N$, $R(0) = 0$, and $D(0) = 0$. In (2.5) and (3.1), the transmission rate, $\beta$, is set to $0.477 = \max_{t\in [0, T]}\zeta(t)$ since $\zeta(t) = \beta(1-\tilde u_1(t))$, where $\tilde u_1(t)$ is a "real-life" non-optimal social distancing control ($0\le \tilde u_1(t) < 1$ for $t\in [0, T]$).

As we look at Figures 2 and 5, it is important to keep in mind that at the start of the study period, in the state of Alabama, the percentage of fully vaccinated people was 33.2%, while in the state of Maryland it was 57%. At the end of the study period, these numbers were 46% and 67.2%, respectively. Hence it comes at no surprise that the daily number of infected people in Alabama (the "reconstructed" scenario) is higher than in Maryland. Figure 5, similar to Figure 2, underscores the importance of disease control in epidemic management. The figure shows a frightening number of infected and deceased people in a hypothetical uncontrolled environment. On the other hand, the "social distancing control" curves in Figure 5 prove the efficiency of this nonmedical form of disease prevention. The "reconstructed" scenario in Figure 5 convincingly shows that the real-life control measures in the state of Maryland worked well.

Figure 6 for Maryland is consistent with what was observed for the state of Alabama in Figure 3. To contain the spread of the virus, social distancing must be strictly enforced at the outset of a new strain. However, this form of control is not sustainable in the long run as the cost begins to take its toll. More than likely, the oscillating behavior of the social distancing control, $u_1(t)$, corresponding to the fourth cost function, $c_{1, 4}$, at the start of the interval is due to unavoidable instability of parameter estimation. As evident from the second graph in Figure 6, similar to the case of Alabama, $u_1(t)$ associated with $c_{1, 4}$ may not be a global minimum. That is why it is so crucial to consider several cost functions to ensure reliable practical recommendations.

Figure 7 illustrates that in the state of Maryland, with optimal social distancing control, the number of new cases and deaths would be significantly reduced during the first 100 days of the Delta strain. However, in the last 40 days, the trend is the exact opposite. Again, this highlights the importance of other control measures, such as vaccination and antiviral treatments, when social distancing inevitably becomes less aggressive over time.

4.   Vaccination control

In this section, we define a mitigation scenario where only vaccination controls are employed, without any social distancing interventions. To accurately capture this scenario, we employ the $S_uS_vI_uI_vRD$ model with control $u_2(t)\in \mathcal{U}$ factored into the vaccination rate. The higher the value of the control, the stricter the vaccination measures, hopefully resulting in a lower number of cases and deaths. Without a social distancing control, the disease transmission rate, $\beta$, is assumed to be at its maximum value, defined as $\beta = \max_{t\in [0, T]}\zeta(t)$. Hence, the $S_uS_vI_uI_vRD$ model with vaccination control, $u_2(t)$, is as follows

In system (4.1), $\nu$ stands for the pre-estimated vaccination capacity, which is set to $1/7$ in our experiments implying that the entire population of the state can potentially be vaccinated in one week once vaccine becomes available for the general population. Of course, in reality, it is impossible to vaccinate everyone, which underlines the importance of the condition $\lim_{u \rightarrow 1^{-}} c_{i, j}(u) = \infty$, guaranteeing that the cost of control becomes prohibitive as $u$ approaches the upper bound of the control range. Recall that for the cost functions $c_{i, j}(u)$, $i = 1, 2$, defined in (2.8), the above assumption is fulfilled when $j = 1, 3, 4$.

In the "vaccine control" scenario, the normalized force of infection that one needs to minimize is equal to $\beta\, \frac{S_u(t)+(1-\alpha)S_v(t)}{1-D(t)}(I_u(t)+I_v(t)),$ leading to the objective functional in the form:

Taking into account Eqs (2.7), (4.1), and (4.2), we arrive at the following Hamiltonian, which is a critical component of the regularized numerical algorithm aimed at estimating the optimal vaccination strategy, $u_2(t)$, for the "vaccine control" problem:

where $\mathbf{p}: = [p_1, p_2, p_3, p_4, p_5, p_6]^\top$. From Pontryagin's Minimum Principle ^[33,83], one concludes that $u_2 = \operatorname*{arg\, min}_{v \in \mathcal{U}} H({\mathbf x}, v, {\mathbf p}),$ subject to state system (4.1) and costate system (2.6):

and $\mathbf{p}(T) = [−1, −1, 0, 0, 0, 0]^\top$. Thus, for the "vaccine control" scenario, one obtains the following $2^{\text{nd}}$-order algorithm for nonlinear constrained minimization:

In Algorithm 2, $F(\theta)$ is a discrete analog of the partial derivative of the Hamiltonian, $H({\mathbf x}, u_2, {\mathbf p})$, with respect to $u_2$, $J(\theta)$ is the Jacobian of $F(\theta)$, $I$ is the identity matrix in the solution space, $\varrho$ is the step size, and $\omega_k$ is the regularization sequence. The derivative of $H({\mathbf x}, u_2, {\mathbf p})$ with respect to $u_2$ exists, since $c_{2, j}(u)$, $j = 1, 2, 3, 4$, are twice continuously differentiable by our assumption. As in the previous section, in all our experiments, shifted Legendre polynomials were used to project the control function, $u_2(t)$, onto a finite dimensional subspace with $\theta$ being a vector of expansion coefficients. MATLAB's built-in function "ode15s" was employed to solve both ODE systems, (4.1) and (4.4), while "lsqnonlin" implemented the Levenberg-Marquardt optimization procedure.

In what follows, we present numerical results comparing the "no control", the "reconstructed", and the "vaccine control" scenarios for managing the COVID-19 Delta variant in the states of Alabama and Maryland. The "reconstructed" (or "real-life") scenario is described by system (2.1), where pre-estimated parameter values for Alabama and Maryland are the same as in Section 3. In the "reconstructed" or "real-life" scenario, we assume that both social distancing and vaccination controls are present, but their implementation is not optimal and mimics real-life interventions put in place from 7/9/2021 to 11/25/2021. The hypothetical "no control" scenario is given by model (2.5), where neither social distancing nor vaccination control is applied. The hypothetical "vaccine control" scenario (4.1) represents the case where optimal vaccination control is implemented in the absence of social distancing measures (see Algorithm 2).

Figure 8 illustrates the disease progression in Alabama from July 9, 2021 to November 25, 2021, highlighting the significant impact of vaccination control throughout the study period. The results suggest that the optimal vaccination strategy involves quickly vaccinating a large portion of the population - over 90%. Numerical experiments demonstrate that vaccinating people early on leads to a substantial decline in the daily number of infected individuals. This translates into a considerable reduction in COVID-related deaths during the same period of time, which further confirms the benefits of early vaccination.

Figure 8 also underscores that the "real-life" preventive measures, including social distancing and vaccination, were very beneficial. At the same time, vaccinating a higher percentage of people before or right after the start of the Delta strain could have prevented even more infections and further depleted the number of deaths. Figure 9 reinforces the message that optimal vaccine control must be applied very early and nearly to its full capacity. Over time, the intensity of vaccine control, understandably, needs to go down due to the lasting benefits of vaccination.

All control strategies shown in Figure 9 adhere to feasible constraints (without these constraints being imposed in the algorithm) for all values of $t\in [0, T]$, with the exception of $u_2(t)$, corresponding to the cost function $c_{2, 2}(u) = u^2$, for which both upper and lower bounds need to be enforced. Thus, Figure 9, once again, demonstrates the importance of the condition $\lim_{u \rightarrow 1^{-}} c_{i, j}(u) = \infty$. The second graph in Figure 9 shows that control strategies associated with $c_{1, 2}$ and $c_{1, 4}$ are not likely to be global minima. Similar to the case of "social distancing control", advantages of the $4^{th}$ control, seen in Figure 8, may be due to insufficient cost reduction.

Examining Figures 4 and 10, one can clearly see a major difference in disease dynamics for "social distancing control" and "vaccine control" scenarios. In the case of optimal "social distancing control" (Figure 4), new COVID cases and deaths are surging toward the end of the study period calling for immediate further interventions. However, the optimal "vaccine control" strategy (Figure 10) reliably "flattens the curve" for the entire interval of time, demonstrating a lasting positive impact of prompt early vaccination.

Figures 11 and 12 make essentially the same case for the optimal "vaccine control" scenario as Figures 8 and 9, that is, for the best outcome, vaccination needs to be done early and it has to include as many people as possible. In terms of the algorithm, Figures 9 and 12 reiterate that condition $\lim_{u \rightarrow 1^{-}} c_{i, j}(u) = \infty$ is critical. For solution $u_2(t)$, corresponding to the cost function $c_{2, 2}(u) = u^2$, to generate a practically relevant control strategy, the bounds of the feasible set need to be enforced in the course of optimization. Once again, the second graph in Figure 12 illustrates that control strategies associated with $c_{1, 2}$ and $c_{1, 4}$ are not likely to be global minima.

The long-term efficiency of the optimal vaccination control is even more pronounced if one compares Figures 7 and 13 for the state of Maryland. In Figure 7, showing the hypothetical social distancing intervention, one can observe a catastrophic increase in daily new cases and deaths in the last 40 days of the study interval. At the same time, for the "reconstructed" scenario in Maryland, in the last 40 days, the epidemic is safely contained. Figure 13, contrasting the optimal "vaccine control" and the "reconstructed" real-life progression, makes a convincing case for the optimal vaccination strategy, which is clearly superior to the "reconstructed" and "social distancing" scenarios by a very large margin.

5.   Both controls concurrently

In this section, we explore a general scenario where both social distancing and vaccination controls are implemented concurrently. To model this optimal control problem, we employ the $S_uS_vI_uI_vRD$ system (2.3) introduced in Section 2. In system (2.3), a control $u_1(t)$ is factored into the transmission rate to account for social distancing and a control $u_2(t)$ is factored into the vaccination rate to account for movement of individuals between susceptible unvaccinated and vaccinated compartments. To balance the pros and cons of optimal controls, $u_1(t)$ and $u_2(t)$, we minimize the objective functional (2.4) subject to system (2.3). As it follows from (2.7), (2.3), and (2.4), the Hamiltonian for the "both controls concurrently" scenario is given by the equation:

Taking into account Pontryagin's Minimum Principle ^[33,83], we recall that the costate system of equations, corresponding to our controlled biological model (2.3), satisfies

This gives rise to the following ODE system, $\textbf{p}(T) = [−1, −1, 0, 0, 0, 0]$:

To minimize (5.1) subject to Eqs (2.3) and (5.2), we propose the Levenberg-Marquardt optimization algorithm:

In Algorithm 3, $F(\theta)$ is a discrete analog of $\partial_{{\mathbf u}}H({\mathbf x},$u$, {\mathbf p})$, $J(\theta)$ is the Jacobian of $F(\theta)$, $I$ is the identity matrix in the solution space, $\varrho$ is the step size, and $\omega_k$ is the regularization sequence. We point out that $\partial_{{\mathbf u}}H({\mathbf x}, {\mathbf u}, {\mathbf p})$ exists, since $c_{i, j}(u)$, $i = 1, 2$, $j = 1, 2, 3, 4$, are twice continuously differentiable by our assumption. In all our experiments, shifted Legendre polynomials were used to project the control function, ${\mathbf u}$, onto a finite dimensional subspace with $\theta$ being a vector of expansion coefficients. Matlab built-in function "ode15s" was employed to solve both ODE systems, (2.3) and (5.2), while "lsqnonlin" implemented the trust-region optimization procedure.

To illustrate the efficiency of the "both controls concurrently" intervention strategy, we present numerical experiments with real data for the SARS-CoV-2 Delta variant of the COVID-19 pandemic in Alabama and Maryland from July 9, 2021 to November 25, 2021 ^[78]. In this section, the "both controls concurrently", "reconstructed", and "no control" scenarios are compared. In the "reconstructed" or "real-life" scenario (2.1) we assume that both social distancing and vaccination controls are present, but their implementation is not optimal and mimics real-life interventions put in place from 7/9/2021 to 11/25/2021. The hypothetical "both controls concurrently" scenario (2.3), represents the case where social distancing and vaccination controls are optimized in the sense of the objective functional (2.4).

The pre-estimated parameters for Alabama in these two cases are $N = 5,031,362$, $\delta_r = 1/90$, $\gamma_{u, r} = (1-0.005)/10$, $\gamma_{u, d} = 0.005/18.5$, $\gamma_{v, r} = (1-0.005/12.7)/10$, $\gamma_{v, d} = 0.005/18.5/12.7$, $\delta_v = 0$, and $\alpha = 0.8$ ^[76,77], while the time-dependent transmission rate, $\zeta(t)$, for (2.1) and case reporting rate, $\psi$, are reconstructed from CDC data ^[78] on daily new infections and deaths by a regularized optimization algorithm. The reconstructed value of $\psi$ is equal to 0.154 (95%CI:[0.149, 0.159]) ^[77]. In (2.1), for the state of Alabama, the average vaccination rate, $\mu$, is equal to $0.0009143$, while the vaccination capacity, $\nu$, in (2.3) is set to $1/7$. The initial values for systems (2.3) and (2.1) are $S_u(0) = 3,402,668/N$, $S_v(0) = 1,626,323/N$, $I_u(0) = 1584/N$, $I_v(0) = 787/N$, $R(0) = 0$, and $D(0) = 0$.

The hypothetical "no control" scenario is given by model (2.5), where neither social distancing nor vaccination control is applied. The initial values in this case are $S_u(0) = (3,402,668+1,626,323)/N$, $I_u(0) = (1584+787)/N$, $R(0) = 0$, $D(0) = 0$. In (2.5) and (3.1), the constant transmission rate, $\beta$, is set to $0.416 = \max_{t\in [0, T]}\zeta(t)$ since $\zeta(t) = \beta(1-\tilde u_1(t))$, where $\tilde u_1(t)$ is a "real-life" non-optimal social distancing control ($0\le \tilde u_1(t) < 1$ for $t\in [0, T]$), as mentioned in Section 2 above.

Figure 14, showing disease progression in the state of Alabama, illustrates superior efficiency of the hypothetical optimal "both controls concurrently" strategy. A comparison of Figures 9 and 15 demonstrates that in the presence of social distancing mitigation measures, the vaccination campaign at the start of the strain can be less rigorous. Figures 10 and 16 underscore that, in combination, optimal social distancing and vaccination controls achieve better results than "vaccine control" without either control being enforced too aggressively. Figures 9 and 15 show that after the initial push, which keeps the epidemic reliably contained, both social distancing and vaccination controls can be scaled down very quickly saving considerable resources. This highlights the importance of simultaneously applying multiple control strategies to combat disease outbreaks and to minimize their negative consequences.

The state-specific parameters for the Maryland "reconstructed" scenario, described by system (2.1), are $N = 6,173,205$ and $\mu = 0.0007286$ ^[76,77]. The time-dependent transmission rate, $\zeta(t)$, and case reporting rate, $\psi$, are estimated from CDC data ^[78] on daily new infections and deaths by a regularized optimization algorithm. The reconstructed value of $\psi$ for Maryland is equal to 0.182 (95%CI:[0.172, 0.192]). In ODE systems (2.1) and (2.3), the initial values for the coordinates of $\mathbf{x}: = [S_u, S_v, I_u, I_v, R, D]^\top$ in the state of Maryland are $S_u(0) = 2,727,503/N$, $S_v(0) = 3,445,221/N$, $I_u(0) = 207/N$, $I_v(0) = 274/N$, $R(0) = 0$, and $D(0) = 0$. ^[77]. As before, the hypothetical "no control" scenario is given by model (2.5) with no control applied. The initial values for the variables of $\mathbf{x}: = [S_u, I_u, R, D]^\top$ are $S_u(0) = (3,402,668+1,626,323)/N$, $I_u(0) = (1584+787)/N$, $R(0) = 0$, and $D(0) = 0$. In (2.3) and (2.5), the constant transmission rate, $\beta$, is set to $0.477 = \max_{t\in [0, T]}\zeta(t)$.

Figures 17–19 emphasize that applying social distancing and vaccine controls in combination is the most powerful way to contain the outbreak. By comparing "both controls concurrently" scenarios for Alabama and Maryland in Figures 15 and 18, respectively, one can notice that in both states, more emphasis is placed on vaccination than on social distancing. However, the two controls quickly decrease at a near-exponential rate due to a sustainable positive effect of vaccination.

As in all previous cases, the $4^{th}$ control may not be a global minimum of the Hamiltonian as suggested by the second graphs in Figures 15 and 18. The first three controls, including the one that corresponds to $c_{1, 2}(u) = c_{2, 2}(u) = u^2$, appear to be global minima, and no bounds need to be enforced by the optimization algorithm.

6.   Conclusions and future plans

Control and prevention strategies are indispensable tools for managing the spread of infectious diseases. This paper examines biological models for the post-vaccination stage of a viral outbreak that integrate two important mitigation tools: social distancing, aimed at reducing the disease transmission rate, and vaccination, which boosts the immune system. Five different scenarios of epidemic progression are considered: (ⅰ) the "no control" scenario, reflecting the natural evolution of a disease without any safety measures in place, (ⅱ) the "reconstructed" scenario, representing real-world data and interventions, (ⅲ) the "social distancing control" scenario covering a broad set of behavioral changes, (ⅳ) the "vaccine control" scenario demonstrating the impact of vaccination on epidemic spread, and (ⅴ) the "both controls concurrently" scenario incorporating social distancing and vaccine controls simultaneously. By investigating these scenarios, we provide a comprehensive analysis of various intervention strategies, offering valuable insights into disease dynamics.

Figures A1–A8 in the Appendix section compare the evolution of state variables, $\hat S_u$, $\hat S_v$, $\hat I_u$, $\hat I_v$, $\hat R$, $\hat D$, for "social distancing control", "vaccine control", and "both controls concurrently" mitigation frameworks in Alabama and Maryland versus the corresponding state variables reconstructed from real data from July 9, 2021 to November 25, 2021 ^[78]. Tables A1–A8 in the Appendix compare the daily number of infected people, $\hat I_u+\hat I_v$, for the same controls and over the same period of time. The main findings of our numerical study are as follows:

● Strict social distancing at the early ascending stage of a new virus wave has immediate positive impact and saves a lot of lives. Early on, optimal social distancing is more efficient than optimal vaccination. However, social distancing is not sustainable in the long run since the cost of this extreme control measure eventually takes its toll. Thus, the intensity of social distancing control begins to slow down after about three months. When social distancing policies are relaxed, new incidence cases and daily new deaths surge. This hypothetical surge is particularly noticeable in Maryland as compared to real data (see Figure 7).

● The positive impact of optimal vaccination is not immediate. But once vaccination takes effect, it works exceptionally well. Experiments show that vaccine control can be safely scaled down after the initial rapid action without any risk of cases increasing or more people dying toward the end of the study window.

● Given the unique individual properties of social distancing and vaccination, the best outcome is achieved when the two controls are applied concurrently. Together, these interventions complement each other. They quickly "flatten the curve" and prevent infections from rebounding in the last 40 days. After prompt early response, both controls quickly decrease at a near-exponential rate due to the sustainable positive effect of vaccination. In all scenarios considered, the results obtained with $c_{i, 4}(u) = - u\ln(1-u)$ should be taken with a grain of salt, since the corresponding control function may not be a global minimum of the Hamiltonian.

● To reconstruct unknown parameters for the optimal control problem, quantification of uncertainty related to noise in the reported data was carried out by refitting model (2.1) to $M = 100$ additional data sets assuming a Poisson error structure. The resulting $M$ best-fit parameter sets were used to find the 95% confidence intervals and to estimate the mean values ^[76,77]. In our future work, we plan to replace system (2.1) with a stochastic compartmental model governed by Wiener processes, which takes into account the uncertainty of the disease transmission, incubation period, and variability of detection.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

A. S. is supported by NSF awards 2011622 and 2409868 (DMS Computational Mathematics).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix