Optimal protection and vaccination against epidemics with reinfection risk

1.   Introduction

As observed during the COVID-19 pandemic, if left unchecked, infectious diseases potentially spread across the entire planet in the span of a few weeks and cause significant damage in terms of mortality and life-long impairments. In addition, the emergence of different variants may lead to reinfection, once the initial immunity weakens over time. In such situations, policy makers impose restrictions on individuals in the form of social distancing and mandatory mask usage. Additionally, they administer vaccines which offer partial immunity against the disease. However, such interventions have a significant social and economic cost, and it is important to strike the right balance among the different options that are available. In this regard, dynamical systems and the optimal control theory have emerged as promising tools that provide policy-makers with appropriate guidelines and insights into mitigating epidemics (see, e.g., ^[1,2]). Additionally, the optimal control of fractional-order systems has been used for spreading processes ^[3].

Starting from seminal works such as ^[4,5], there have been numerous investigations on optimal control of epidemiological processes, which largely consider compartmental dynamical models of epidemic evolution ^[1]; see ^[2] for a recent review. A majority of the past efforts have been directed towards optimal protection in the context of Susceptible-Infected-Recovered (SIR) epidemics and its variants (see, e.g., ^[6,7,8,9]). More recent papers ^{[10,11,12,13]} have considered vaccinations as additional control input (in addition to protection or social distancing measures). These works assume that the running cost is quadratic in the control input. However, it is natural to assume that the cost (of vaccination or protection) is proportional to the magnitude of the control input or the fraction of the population on which the input is administered. A few additional papers (see, e.g., ^[14,15]) have investigated the use of optimal control techniques, when the population size is dynamically changing. Other related approaches are also explored in ^[16,17,18]. The methodology for containing COVID-19 Delta strain spread was explored in ^[19]. The authors included asymptomatic agents and captured the notion of an imperfect vaccination in their model. It is well established that fractional order optimal controls have advantages in the form of greater flexibility and higher accuracy over the classical integer order controls. The fractional order models of COVID-19 and other diseases were thoroughly explored in ^[20,21,22].

There have been limited investigations into epidemiological models where recovery does not give permanent immunity, and hence, reinfection is also possible as a result. In addition, even past works that assumed a cost functional that was linear in the control input which led to a bang-bang optimal control structure, the possibility of the existence of singular arcs and singular control inputs is not often examined in a rigorous manner (the work ^[23] is a notable exception in this regard). Nevertheless, in practice, it is important to characterize the possibility of singular control inputs in order to provide insights into policy-making decisions, thus informing the authorities of the expected impact of either imposing or relaxing interventions.

The motivation for this work is to establish the existence of non-singular optimal policies to control the spread of epidemics via limited vaccination and protective measures by solving an optimal control problem considering a running cost that is linear with respect to the input variables in an epidemic model with reinfection risk. Our setting differs from most prior studies on the optimal control of epidemics that assume the objective function to be superlinear in the control inputs, which leads to a simpler analysis, and the issue of singular inputs can be avoided. For example, in ^[24], the authors used the SIR model where the objective function was quadratic in the control inputs. However, it is more reasonable to consider the running cost to be linear with respect to the input variables; indeed, the cost of vaccination (and other protection measures) is directly proportional to the fraction of the population being vaccinated (or adopting protective measures). While some studies, such as ^[25], assumed running costs that were linear in the control inputs, they focused on bang-bang controls without ruling out the possibility of singular controls. In this work, we consider a generalized epidemiological model that incorporates both recovery and reinfection (similar to observations made during COVID-19), specifically the susceptible-infected-recovered-infected (SIRI) epidemic model (see e.g., ^[26]). Our model includes both non-pharmaceutical and medical resources as inputs, and the running cost is assumed to be linear in these control variables. Additionally, during COVID-19, we observed higher reinfection rates due to variants such as Delta and Omicron, which supports the focus on compromised immunity in this work. Under appropriate assumptions, we specifically exclude the possibility of simultaneous singularities and analytically prove that vaccination control does not exhibit singularities.

In the SIRI model, the rate of reinfection is different from the rate of initial infection, with higher values indicating compromised immunity and a smaller rate of reinfection indicating a partial immunity imparted by the disease and/or vaccinations. As analyzed in ^[27], it was assumed that vaccinations were only available for the susceptible sub-population who transit to the recovered compartment, thus reflecting the fact that vaccinations impart a certain degree of protection for the short term, but not complete immunity (a similar phenomenon was also observed during the COVID-19 pandemic).

The main contributions of the paper are as follows. We analyze the optimality conditions for the associated optimal control problem and rule out the existence of a simultaneous singularity of both control inputs on the SIRI epidemiological model in the scenario of compromised immunity. Then, we carry out a detailed investigation regarding the singularity of the vaccination control input, and under sufficient conditions, we show that it does not admit a singular arc (i.e., the vaccination-optimal control is always at one of two possible extreme admissible values). A theoretical analysis that provides valuable insight on the vaccination-control input being non-singular (also known as bang-bang or, on-off control) is essential, since bang-bang control is often considered a more appropriate intervention in practical epidemiological and clinical settings (see, e.g., ^[28,29]). Additionally, we demonstrate epidemic evolution under optimal control inputs for a numerical case study and show the relative impact of vaccination and protection in an epidemic containment.

The remainder of the paper is organized as follows: the controlled SIRI epidemiological model is introduced in Section 2, where we also prove the existence and local asymptotic stability of its equilibrium points; the optimal control problem is presented in Section 3 and the structural properties of the optimal control inputs are established; the non-existence of singular arcs in the structure of the candidate optimal control corresponding to vaccination is established in Section 4; numerical results depicting the evolution of the epidemic under the optimal control inputs are presented in Section 5; and we conclude in Section 6 with a discussion on possible directions for future research.

2.   Controlled SIRI epidemiological model

Motivated by the COVID-19 pandemic, we consider the SIRI epidemiological model, which was introduced in ^[26]. In this setting, an individual remains in one of three possible states: susceptible $(\mathtt{S})$, infected $(\mathtt{I})$, or recovered $(\mathtt{R})$. However, recovery is not permanent, and recovered individuals also become potentially infected again upon contact with infected individuals. The rate at which a susceptible (respectively, recovered) individual becomes infected upon contact with infected individuals is denoted by $\beta > 0$ (respectively, $\hat{\beta} > 0$). In general, $\beta$ is assumed to be different from $\hat{\beta}$. When $\hat{\beta} < \beta$, the reinfection rate is smaller than the rate of new infection, which indicates that recovery imparts a partial immunity against future infection. Similarly, $\hat{\beta} > \beta$ indicates a compromised immunity following the initial infection. Finally, $\gamma > 0$ represents the rate at which an infected individual recovers, which is referred to as recovery rate. The various state transitions are depicted in Figure 1.

We consider two types of control inputs (which are assumed to be essentially bounded Lebesgue measurable functions): $u_{\texttt{V}}$, which captures the rate at which susceptible individuals are vaccinated, and $u_{\mathtt{P}}$, where $1 - u_{\mathtt{P}}$ captures the effective rate of social distancing or protection adoption by individuals in the disease states $\texttt{S}$ and $\texttt{R}$. As a consequence of the above control inputs, the resulting controlled SIRI epidemic dynamical equations are given by the following:

where the state variables $x_{\texttt{S}}(t) \in [0, 1], x_{\texttt{I}}(t) \in [0, 1]$ and $x_{\texttt{R}}(t) \in [0, 1]$ denote the instantaneous proportion of individuals in each of the three epidemic states $\mathtt{S}, \mathtt{I}$, and $\mathtt{R}$. Henceforth, unless required, we suppress the dependency of the states and control inputs on time $t$ for an improved readability.

Remark 2.1. The biological significance to controlling the epidemic spread is that our analysis accounts for reinfection in individuals, which aligns with the characteristics of several infectious diseases, such as COVID-19, that only confers a short-term immunity. In addition, during the COVID-19 pandemic, it was demonstrated that the reinfection rates, particularly due to variants such as Delta and Omicron, exceed the initial infection rates ^[30]. Motivated by this observation, we later assume that $\hat{\beta} > \beta$, which implies that getting infected compromised the immunity. These characteristics are not captured by the classical epidemic models, such as the SIR model.

Remark 2.2. Note that the term $\beta x_{\mathtt{S}} x_{\mathtt{I}} u_{\mathtt{P}}$ represents the fraction of the susceptible population who do not adopt any protection and get infected, while the term $x_{\mathtt{S}} u_{\mathtt{V}}$ represents the fraction of the susceptible population who opt for vaccination and move to the recovered state. Note that such individuals may become infected again in future (i.e., vaccination does not impart a permanent immunity against future infection). Similarly, the term $\hat{\beta} x_{\texttt{R}} x_{\texttt{I}} u_{\mathtt{P}}$ captures the fraction of the recovered population who do not adopt any protection and get reinfected, and the term $\gamma x_{\texttt{I}}$ is the fraction of the infected population who naturally recover. The above dynamics satisfies $\dot{x}_{\texttt{S}}(t) + \dot{x}_{\texttt{I}}(t) + \dot{x}_{\texttt{R}}(t) = 0$ for almost every instant of time $t$, and since the states represent fractions of the population, they also satisfy $x_{\texttt{S}}(t) + x_{\texttt{I}}(t) + x_{\texttt{R}}(t) = 1$ for every instant of time $t$, when the initial state vector also satisfies this condition (for details, see Lemma 2.1).

In our model, the control inputs include behavioral measures (represented by $1 - u_{\mathtt{P}}$), such as protective behaviors or social distancing, and medical interventions (represented by $u_{\mathtt{V}}$), such as vaccination. Thus, our model accounts for both medical and non-medical interventions available during an epidemic. We impose limitations on both types of inputs to prevent trivial solutions that might arise from an unlimited supply of protection and vaccination. In the above setting, $u_{\mathtt{P}} = 0$ implies that the susceptible or recovered individuals adopt complete a protection and they do not bear the risk of getting infected. In order to rule out this impractical corner case, we assume that $u_{\mathtt{P}}$ is always bounded from below by a lower bound $u_{\texttt{Pmin}} > 0$. Additionally, we assume that $u_{\texttt{P}} \leq 1$ with the upper bound chosen to signify that $\beta$ and $\hat{\beta}$ denote the infection rates in the absence of any protective action. In addition, we assume that the vaccination rate satisfies $0 \leq u_{\texttt{V}} \leq u_{\texttt{Vmax}} < 1$, where we have limited the upper threshold by excluding $1$, as $u_{\mathtt{V}} = u_{\texttt{Vmax}} = 1$ would imply vaccinating the entire susceptible fraction of population in one go, which is not practical.

Remark 2.3. When $\hat{\beta} = 0$ (i.e., recovered individuals do not get reinfected), then the model reduces to the SIR epidemiological model (see ^[31]). Similarly, as mentioned in ^[31], when $\beta = \hat{\beta}$ (i.e., the infection rate of susceptible and recovered individuals coincide), then we recover the Susceptible-Infected-Susceptible (SIS) epidemiological model. Thus, the SIRI epidemiological model studied in this paper is a strict generalization of both the SIS and SIR epidemiological models (see, e.g., ^[26]).

Before stating the optimal control problem studied in this paper, we first establish certain theoretical properties of the controlled SIRI epidemiological model when the control inputs are exogenous constants. First, we investigate the equilibria of its dynamics and their associated stability properties. When $u_{\texttt{V}} = 0$, the dynamics in (2.1) is an instance of the SIRI model without any explicit control input, with the infection rates effectively being $\beta u_{\mathtt{P}}$ and $\hat{\beta} u_{\mathtt{P}}$, respectively. The equilibria and their stability properties follow from analogous results established for the classical SIRI epidemiological model in ^[19]. Therefore, we focus on the case where the constant steady-state inputs are defined by $u_{\texttt{V}} = u_{\mathtt{V}}^{\texttt{eq}}$, where $0 < u_{\mathtt{V}}^{\texttt{eq}} \leq u_{\texttt{Vmax}}$, and $u_{\mathtt{P}} = u_{\mathtt{P}}^{\texttt{eq}}$, where $u_{\texttt{Pmin}} \leq u_{\mathtt{P}}^{\texttt{eq}} \leq 1$.

By equating the right hand side of (2.1) to zero, we observe the existence of two equilibrium points:

(i) The disease free equilibrium point $E_{\mathtt{DFE}}$ for $(x^{\texttt{eq}}_\texttt{S} = 0, x^{\texttt{eq}}_\texttt{I} = 0, x^{\texttt{eq}}_\texttt{R} = 1)$, which always exists;

(ii) The endemic equilibrium point $E_{\mathtt{EE}}$ for $(x^{\texttt{eq}}_\texttt{S} = 0, x^{\texttt{eq}}_\texttt{I} = 1 - \frac{\gamma}{\hat{\beta} u^{\texttt{eq}}_\texttt{P}}, x^{\texttt{eq}}_\texttt{R} = \frac{\gamma}{\hat{\beta} u^{\texttt{eq}}_\texttt{P}})$, which exists when $\gamma < \hat{\beta} u^{\texttt{eq}}_\texttt{P}$.

The following result establishes their local stability properties.

Proposition 2.1 (Local asymptotic stability of the equilibrium points). For the controlled SIRI epidemiological model (2.1) with $u^{ {{\texttt{eq}}}}_ {{\texttt{V}}} > 0$, we have the following:

(i) The disease free equilibrium point is locally asymptotically stable when $\gamma > \hat{\beta} {u^{ {{\texttt{eq}}}}_ {{\texttt{P}}}}$;

(ii) The endemic equilibrium point is locally asymptotically stable when $\gamma < \hat{\beta} {u^{ {{\texttt{eq}}}}_ {{\texttt{P}}}}$.

Proof. Since $x_{\texttt{S}}(t) + x_{\texttt{I}}(t) + x_{\texttt{R}}(t) = 1$ for every instant of time $t$ (see Lemma 2.1), we equivalently consider the dynamics involving only the two state variables $x_\texttt{S}$ and $x_\texttt{I}$, by expressing $x_\texttt{R} = 1 - x_\texttt{S} - x_\texttt{I}$. Thus, the dynamics (2.1) reduces to the following:

the Jacobian matrix of which is given by the following:

First, we investigate the local asymptotic stability of the disease-free equilibrium point. In this case, the Jacobian matrix is given by the following:

Since $u^{\texttt{eq}}_\texttt{V} > 0$, both the eigenvalues of the above matrix, are strictly negative when $\gamma > \hat{\beta} u_{\mathtt{P}}^{\texttt{eq}}$.^* Next, we investigate the local asymptotic stability of the endemic equilibrium point. The Jacobian matrix in this case, is given by

^*Intuitively, when $\gamma < \hat{\beta} u^{\texttt{eq}}_\texttt{P}$, the recovery rate of the infected fraction of population is less than the rate of reinfection of the recovered fraction of population times the fraction of recovered people not adopting protection, leading to the disease free equilibrium being unstable.

It is easy to see that both the eigenvalues of the above matrix, are strictly negative when $\gamma < \hat{\beta} u_{\mathtt{P}}^{\texttt{eq}}$. This concludes the proof. □

The following lemma supports the consideration of the SIRI epidemiological model.

Lemma 2.1 (Positive invariant set of the controlled SIRI epidemiological model). The set $S : = \{(x_{ {{\texttt{S}}}}, x_{ {{\texttt{I}}}}, x_{ {{\texttt{R}}}}) : (x_{ {{\texttt{S}}}}, x_{ {{\texttt{I}}}}, x_{ {{\texttt{R}}}}) \in [0, 1] \times [0, 1] \times [0, 1]\}$ is positively invariant, with respect to any unique global solution for the SIRI epidemic dynamics (2.1).

Proof. First, we show that local solutions exist and are also unique on a sufficiently small time-interval for the SIRI epidemic dynamics (2.1). To this end, fix a sufficiently small real number $\varepsilon > 0$, and for any given Lebesgue measurable control inputs $\textbf{u} = (u_{\mathtt{P}}, u_{\mathtt{V}}): [0, \varepsilon]\rightarrow [u_{{\texttt{Pmin}}}, 1]\times [0, u_{{\texttt{Vmax}}}]$, let us rewrite (2.1) as follows:

where $\textbf{x} : = (x_{\mathtt{S}}, x_{\mathtt{I}}, x_{\mathtt{R}}) \in \mathbb{R}^{3}$ is the state vector and the function $\textbf{G}: \mathbb{R} \times \mathbb{R}^3 \rightarrow \mathbb{R}^3$ is given by the following:

Now, let $K : = \{(t, \textbf{x}) : 0 \leq t \leq \varepsilon, |\textbf{x} - \textbf{x}_{0}| \leq \varepsilon\}$ be a cylinder in $\mathbb{R} \times \mathbb{R}^{3}$.^† It is now easy to verify that the function $\textbf{G}$, satisfies the following conditions:

^†The norm of a vector $\textbf{x} \in \mathbb{R}^{3}$, is denoted by $|\textbf{x}|$.

(i) For almost every $t \in [0, \varepsilon]$, the mapping $\textbf{x}\mapsto \textbf{G}(t, \textbf{x})$ is continuous, and for every $\textbf{x} \in \mathrm{B}_{\varepsilon}(\textbf{x}_{0})$, the mapping $t\mapsto \textbf{G}(t, \textbf{x})$ is Lebesgue measurable; ^‡

^‡The closed ball of radius $r > 0$ in $\mathbb{R}^{3}$, centered at $\textbf{x} \in \mathbb{R}^{3}$, is denoted by $\mathrm{B}_{r}(\textbf{x})$.

(ii) There exists a constant $C_{K} > 0$ such that:

holds for almost every $t \in [0, \varepsilon]$ and for every $\textbf{x} \in \mathrm{B}_{\varepsilon}(\textbf{x}_{0})$. Moreover, there also exists a constant $L_{K} > 0$ such that the following inequality:

holds for almost every $t \in [0, \varepsilon]$ and for every $\textbf{x}, \textbf{y} \in \mathrm{B}_{\varepsilon}(\textbf{x}_{0})$.

Indeed, to verify the first claim in item (ii) stated above, one can obtain the following:

which holds for almost every $t \in [0, \varepsilon]$ and for every $\textbf{x} \in \mathrm{B}_{\varepsilon}(\textbf{x}_{0})$. Keeping in mind, the fact that the functions $\textbf{f}, \textbf{g}_{1}, \textbf{g}_{2}$ are of class $C^{1}$, one can now obtain the desired result by invoking Weierstrass' theorem. To verify the second claim in item (ii) stated above, one can obtain the following inequality:

which holds for almost every $t \in [0, \varepsilon]$ and for every $\textbf{x}, \textbf{y} \in \mathrm{B}_{\varepsilon}(\textbf{x}_{0})$. Keeping in mind, the facts that the functions $\textbf{f}, \textbf{g}_{1}, \textbf{g}_{2}$ are of class $C^{1}$ and also that the set $\mathrm{B}_{\varepsilon}(\textbf{x}_{0})$ is compact and convex, one can now obtain the desired result. By appropriately modifying some of the steps given in the proof of [32, Theorem 2.2.1], one can now deduce that local solutions exist for (2.1) on the time-interval $[0, \epsilon]$. In addition, from [32, Theorem 2.1.3], one can also deduce the uniqueness of such local solutions of (2.1) on the time-interval $[0, \hat \epsilon]$, where $0 < \hat \varepsilon \leq \varepsilon$.

Next, we verify the positive invariance of the set $S$ with respect to the unique local solution $\textbf{x} = (x_{\mathtt{S}}, x_{\mathtt{I}}, x_{\mathtt{R}}): [0, \hat \varepsilon] \rightarrow \mathbb{R}^3$ of the SIRI epidemic dynamics (2.1). To this end, from (2.1), we have that

for all $t \in [0, \hat \varepsilon]$. Now, observe that the initial conditions $x_{\mathtt{S}}(0), x_{\mathtt{I}}(0), \text{ and } x_{\mathtt{R}}(0)$ represent the initial fractions of the population who are susceptible, infected, and recovered, respectively, and therefore satisfy the following constraints: $x_{\mathtt{S}}(0), x_{\mathtt{I}}(0), x_{\mathtt{R}}(0) \geq 0$ and $x_{\mathtt{S}}(0) + x_{\mathtt{I}}(0) + x_{\mathtt{R}}(0) = 1$. Since, $x_{\mathtt{S}}(0), x_{\mathtt{I}}(0)\geq 0$, it follows from the first two equations in (2.2) that $x_{\mathtt{S}}(t), x_{\mathtt{I}}(t)\geq 0$ for all $t \in [0, \hat \varepsilon]$, and since $x_{\mathtt{R}}(0) \geq 0$, $u_{\mathtt{V}}(t)\geq 0$ for all $t \in [0, \hat \varepsilon]$ and $\gamma \geq 0$, it now follows from the third equation in (2.2) that $x_{\mathtt{R}}(t)\geq 0$ for all $t \in [0, \hat \varepsilon]$. Now, from the dynamics (2.1), it follows that the relation $\dot{x}_{\mathtt{S}}(t) + \dot{x}_{\mathtt{I}}(t) + \dot{x}_{\mathtt{R}}(t) = 0$ is satisfied for almost every $t \in [0, \hat \varepsilon]$, which in turn implies that the relation $x_{\mathtt{S}}(t) + x_{\mathtt{I}}(t) + x_{\mathtt{R}}(t) = 1$ is satisfied for all $t \in [0, \hat \varepsilon]$. Overall, the unique local solution $\textbf{x} = (x_{\mathtt{S}}, x_{\mathtt{I}}, x_{\mathtt{R}}): [0, \hat \varepsilon] \rightarrow \mathbb{R}^3$ of the SIRI epidemic dynamics (2.1), satisfies the following constraints: $x_{\mathtt{S}}(t), x_{\mathtt{I}}(t), x_{\mathtt{R}}(t) \geq 0$ and $x_{\mathtt{S}}(t) + x_{\mathtt{I}}(t) + x_{\mathtt{R}}(t) = 1$ for all $t \in [0, \hat \varepsilon]$, from which it follows that $(x_{\mathtt{S}}(t), x_{\mathtt{I}}(t), x_{\mathtt{R}}(t)) \in [0, 1] \times [0, 1] \times [0, 1]$ for all $t \in [0, \hat \varepsilon]$.

Finally, we show that any unique right-maximal solution $\textbf{x} = (x_{\mathtt{S}}, x_{\mathtt{I}}, x_{\mathtt{R}}): [0, T) \rightarrow \mathbb{R}^3$ of the SIRI epidemic dynamics (2.1), where the time $T > 0$, can be globally extended (i.e., it is possible to show that this holds for $T = \infty$). To this end, let us assume that $T < \infty$; then, by appropriately modifying some of the steps given in the proof of [32, Theorem 2.1.4], one can deduce the following relation:

where the set $\Omega : = [0, T) \times \mathbb{R}^{3}$ is an open set in $[0, \infty] \times \mathbb{R}^{3} \subset \mathbb{R} \times \mathbb{R}^{3}$, the notation $\partial \Omega$ denotes its boundary (in the set $[0, \infty] \times \mathbb{R}^{3}$), and the distance of a pair of points $(t, \textbf{y}) \in [0, \infty] \times \mathbb{R}^{3}$ to a set $K \subseteq [0, \infty] \times \mathbb{R}^{3}$ is given by $d((t, \textbf{y}), K) : = \inf_{(\bar t, \bar {\textbf{y}}) \in K} |(t, \textbf{y}) - (\bar t, \bar {\textbf{y}})|$. Using the facts that $S$ is a compact set in $\mathbb{R}^{3}$ and is positively invariant with respect to the unique right-maximal solution $\textbf{x} = (x_{\mathtt{S}}, x_{\mathtt{I}}, x_{\mathtt{R}}): [0, T) \rightarrow \mathbb{R}^3$ of the SIRI epidemic dynamics (2.1), together with the fact that the boundary $\partial \Omega = (\{T\} \times \mathbb{R}^{3}) \cup ([0, T] \times \emptyset)$, we now arrive at a contradiction in view of (2.3). This completes the proof. □

Remark 2.4. Note that the positive invariance of the set $S$ in the proof of Lemma 2.1 can also be shown using Nagumo's theorem adapted to control systems (see, e.g., [33, Theorem $4.11$]).

3.   Optimal control problem

Now, we consider the following optimal control problem:

where $u_{{\texttt{Pmin}}} > 0$ represents the minimum fraction of the susceptible or recovered sub-populations who remain unprotected, and $u_{{\texttt{Vmax}}} < 1$ denotes an upper bound on the fraction of the total population that can be vaccinated for a given time period. The individual weighing terms in the running cost (3.1) are as follows: $c_{\mathtt{P}}$ captures the cost incurred due to the protection adopted by the susceptible and recovered individuals, $c_{\mathtt{V}}$ captures the cost incurred due to vaccination by the susceptible individuals, and $c_{\mathtt{I}}$ is the disease cost or the cost incurred on being infected.

3.1.   Equivalent formulation in Mayer form

In order to exploit the results from the optimal control theory in our subsequent analysis, we convert the optimal control problem defined by (3.1) into the Mayer form. In the Mayer form, the cost functional solely consists of the terminal cost. This requires appending an additional state $x_\texttt{C}$ to our dynamics which captures the running cost, and whose time-evolution satisfies the following dynamics:

for almost every time instant $t$. In addition, we note that any one of the three epidemic states can be expressed in terms of the other two states since the SIRI dynamics (2.1) satisfies $x_\texttt{S}(t)+x_\texttt{I}(t)+x_\texttt{R}(t) = 1$ for every instant of time $t$ and $\dot{x}_\texttt{S}(t)+\dot{x}_\texttt{I}(t)+\dot{x}_\texttt{R}(t) = 0$ for almost every instant of time $t$. As a result, the epidemic dynamics can be expressed in terms of only two state variables. We express $x_\texttt{I} = 1-x_\texttt{S}-x_\texttt{R}$ and omit the variable $x_\texttt{I}$ from the epidemic dynamics. By introducing the state vector $\textbf{z} = ({x}_\texttt{C}, x_\texttt{S}, x_\texttt{R}) \in \mathbb{R}^{3}$, the optimal control problem defined by (3.1) can now be written in the Mayer form as follows:

where the drift and control vector fields are given by the following:

3.2.   Existence of an optimal control input

Now, we establish the existence of a solution for the optimal control problem defined by (3.3). To this end, we leverage Filippov's theorem, which is stated below for the reader's convenience.

Theorem 3.1. (Filippov's theorem, [35, Section 4.5]) Consider a controlled dynamical system:

where $\textbf{x}(t) \in \mathbb{R}^{n}$ and $\textbf{u}(t) \in U \subset \mathbb{R}^{m}$. Assume that the solutions of (3.4) exist on a given time-interval $[0, T]$ for all control inputs $\textbf{u}(\cdot)$. In addition, assume that for every pair $(t, \textbf{x}) \in [0, T] \times \mathbb{R}^{n}$, the set $\{ f(t, \textbf{x}, \textbf{u}) : \textbf{u} \in U\}$ is compact and convex. Then, the reachable set $R^t(\textbf{x}_0)$ is compact for each $t \in [0, T]$.^§

^§The reachable set at time $t > 0$, starting from $\textbf{x}_0 \in \mathbb{R}^{n}$, is defined as follows:

Now, we state the following theorem.

Theorem 3.2 (Existence of an optimal control). There exists a solution for the optimal control problem defined by (3.3).

Proof. In view of Lemma 2.1, it is clear that there exists a unique solution (with respect to any given initial condition and admissible control inputs) defined over the time-interval $[0, T]$ for the dynamics given in the optimal control problem defined by (3.3). Moreover, it is easy to verify that for every $\textbf{z} \in \mathbb{R}^{3}$, the set $\{\textbf{f}(\textbf{z})+\textbf{g}_{\texttt{P}}(\textbf{z}) u_{{\texttt{P}}} +\textbf{g}_{\texttt{V}}(\textbf{z}) u_{{\texttt{V}}} : (u_{\mathtt{P}}, u_{\mathtt{V}}) \in [u_{{\texttt{Pmin}}}, 1] \times [0, u_{\texttt{Vmax}}]\}$ is a compact and convex set in $\mathbb{R}^{3}$. Now, it follows from Theorem 3.1 that the reachable set at time $T$, starting from any given initial condition, is a compact set in $\mathbb{R}^{3}$. Now, the proof can be concluded by invoking Weierstrass' extreme value theorem. □

3.3.   Structure of optimal control inputs

Now, we establish the structure of the candidate optimal control inputs. To this end, we first use Pontryagin's maximum principle to single out the optimal control inputs. The Hamiltonian function which corresponds to the optimal control problem defined by (3.3) is given by the following:

where $\langle \cdot, \cdot \rangle$ denotes the standard inner-product of two vectors in $\mathbb{R}^3$ and, $\pmb\lambda = (\lambda_\texttt{C}, \lambda_\texttt{S}, \lambda_\texttt{R}) \in \mathbb{R}^{3}$ denotes the co-state vector.^¶ For almost every time $t \in [0, T]$, the minimizing control inputs are given by the following:

^¶The absence of the case of the abnormal multiplier being equal to zero for the optimal control problem defined by (3.3) directly follows as a consequence of [34, Corollary 22.3].

where the superscript $(\cdot)^*$ denotes the optimal trajectories, and the co-state dynamics are given by

which satisfies the following terminal boundary condition:

From (3.6), we obtain the following co-state dynamics:

Since the Hamiltonian function is affine with respect to the control inputs, the structure of the optimal control inputs will be governed by the so-called switching functions given by the following (see, e.g., [35, Section 4.4.3]):

The structure of the optimal control inputs $u_{\mathtt{P}}^*$ and $u_{\mathtt{V}}^*$ are given as follows:

where $\star$ denotes the unknown candidate optimal control input, which is also referred to as a singular control input.

When $\phi_{\texttt{Q}} \not \equiv 0$ (i.e., $\phi_{\texttt{Q}}$ is not identically zero on an open time-interval of $[0, T]\subset \mathbb{R}$) for $\texttt{Q} \in \{\texttt{P}, \texttt{V}\}$, then the optimal control inputs $u_{\mathtt{V}}^*$ and $u_{\mathtt{P}}^*$ switch between their respective minimum and maximum admissible values, depending upon the sign of ${\phi}_{\texttt{Q}}$. The control inputs with this property are called bang-bang control inputs. However, we may also encounter a situation in which $\phi_{\texttt{Q}} \equiv 0$ is accompanied by the higher derivatives of $\phi_{\texttt{Q}}$ also vanishing on an open time-interval of $[0, T]\subset \mathbb{R}$, i.e., $\dot{\phi}_{\texttt{Q}} \equiv 0, \; \ddot{\phi}_{\texttt{Q}} \equiv 0, \; \dddot{\phi}_{\texttt{Q}} \equiv 0$, and so on. The control inputs which exhibit such a phenomenon are called singular control inputs (see, e.g., ^[36]).

4.   Non-existence of singular control inputs

An important mathematical tool required for the analysis of a singularity of an optimal control input is the Lie bracket. Let $\textbf{f}$ and $\textbf{g}$ be two continuously differentiable vector fields defined in $\mathbb{R}^n$. Then, for any given $\textbf{x} \in \mathbb{R}^n$, their Lie bracket is defined as follows:

where $D\textbf{f}(\textbf{x})$ and $D\textbf{g}(\textbf{x})$ denote the Jacobian matrices of the vector fields $\textbf{f}$ and $\textbf{g}$, respectively, evaluated at the point $\textbf{x} \in \mathbb{R}^n$.

4.1.   Simultaneous non-singularity of the optimal control inputs $u_{\mathtt{P}}^*$ and $u_{\mathtt{V}}^*$

The existence of bang-bang control inputs (or equivalently non-existence of singular control inputs) is determined by the switching functions and their higher time-derivatives. As previously discussed, a singularity of an optimal control input arises when the switching function, $\phi_{\mathtt{Q}}$ for $\texttt{Q} \in \{\texttt{P}, \texttt{V}\}$ vanishes identically over some time-interval, which is open in $[0, T]$.

First, we examine the possibility of both candidate optimal control inputs being simultaneously singular. It can be shown (see, e.g., [35, Section $4.4.3$]) that the switching functions $\phi_{\mathtt{P}}$ and $\phi_{\mathtt{V}}$ given by (3.9) and (3.10), respectively, have higher order time-derivatives given by the following:

for $\mathtt{Q} \in \{\mathtt{P}, \mathtt{V}\}$. Before we state our result, we introduce the following assumption.

Assumption 4.1. We proceed with the analysis on the class of diseases in which the reinfection rate exceeds the first time infection rate (i.e., $\hat{\beta} > \beta$).

The generalized SIRI epidemiological model discussed in this work takes reinfection into account. The finite-horizon optimal control problem presented in Eq (3.3) involves several model parameters ($\beta, \hat{\beta}, \gamma$), costs ($c_{\mathtt{P}}, c_{\mathtt{I}}, c_{\mathtt{V}}$), and two control inputs ($u_{\mathtt{P}}, u_{\mathtt{V}}$) with defined lower and upper thresholds. Due to the large number of variables and parameters, it is quite challenging to derive a complete characterization for such a model. As a result, we focus on characterizing the control input behaviors for the case that the reinfection rate exceeds the initial infection rate satisfying $\hat{\beta} > \beta$. As mentioned above, during the recent COVID-19 pandemic, it was observed that the reinfection rates, particularly those associated with variants such as Delta and Omicron, were higher than the initial infection rates. The results of this paper are applicable to such immuno-compromising infectious diseases.

Now, we state the following proposition.

Proposition 4.1. Suppose Assumption 4.1 holds; then, the candidate optimal control inputs $u_{\mathtt{V}}^*$ and $u_{\mathtt{P}}^*$ cannot be simultaneously singular on any open time-interval $I \subset [0, T]$.

Proof. As discussed in Section 3.3, a control input exhibits a singularity when the switching function associated with it and its higher derivatives are all identically zero over an open time-interval. In our setting which is comprised of two control inputs, the necessary condition for the existence of a simultaneous singularity of the inputs $u_{\mathtt{V}}^*$ and $u_{\mathtt{P}}^*$ on $I$ requires the following:

to hold for every/almost every $t \in I$, which implies the following:

where we obtain (4.5) from (4.2) and (4.3), with

It follows that (4.5) holds if either of the two conditions is true: the co-state vector identically vanishes (i.e., $\pmb\lambda^*(t) \equiv \textbf{0}$) or the vectors $\textbf{g}_{\mathtt{V}}(\textbf{z}^*(t))$, $\textbf{g}_{\mathtt{P}}(\textbf{z}^*(t))$ and $[\textbf{f} + \textbf{g}_{\mathtt{P}} u_{\mathtt{P}}^*, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t))$ are linearly dependent over $I$. Now, from (3.7) and (3.8), we conclude that $\lambda^*_{\texttt{C}}(t) \equiv 1$ is satisfied for every $t \in I$, which implies that the co-state vector $\pmb\lambda^*(t) \not \equiv \textbf{0}$. Thus, the vectors $\textbf{g}_{\mathtt{V}}(\textbf{z}^*(t))$, $\textbf{g}_{\mathtt{P}}(\textbf{z}^*(t))$, and $[\textbf{f} + \textbf{g}_{\mathtt{P}} u_{\mathtt{P}}^*, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t))$ must be linearly dependent over $I$ for a simultaneous singularity of the inputs to exist. By computing the determinant of the matrix formed by these three vectors, we obtain the following (we have suppressed the explicit dependency of the states and control inputs on time for the sake of brevity):

Suppose the three vectors are linearly dependent. Observe from (2.2) that $x_{\mathtt{S}}^{*}$ is strictly positive for a given bounded time-interval. Similarly, $x_{\mathtt{I}}^{*} = 1 - x_{\mathtt{S}}^{*} - x_{\mathtt{R}}^{*}$ is also non-zero in the endemic case. In addition, $0 < u_{\mathtt{Pmin}} \leq u_{\mathtt{P}}^{*} \leq 1$ implies that the input $u_{\mathtt{P}}^{*}$ is strictly positive. Thus, $x_{\mathtt{S}}^{*{^2}} (1 - x_{\mathtt{S}}^* - x_{\mathtt{R}}^*) u_{\mathtt{P}}^*$ is clearly non-zero on $I$. Now, setting the determinant in (4.7) equal to $0$ and using the relation $x_{\mathtt{S}}^*(t) + x_{\mathtt{R}}^*(t) \neq 0$ results in the following:

Observe from (2.2) that $x_{\mathtt{S}}^{*}(t)$ is an exponentially decreasing function, which remains strictly positive for a given bounded time-interval. Similarly, Lemma 2.1 ensures that $x_{\mathtt{R}}^{*}(t) \geq 0$ holds. The left-hand side of (4.8) corresponds to the state-variable $x_{\mathtt{I}}^*(t)$, which resides in the set $[0, 1]$ by Lemma 2.1. Whereas, the right-hand side is a negative constant under a compromised immunity (i.e., $\hat{\beta} > \beta$), thus implying that the equality (4.8) can never hold. Hence, the three vectors $\textbf{g}_{\mathtt{V}}(\textbf{z}^*(t))$, $\textbf{g}_{\mathtt{P}}(\textbf{z}^*(t))$, and $[\textbf{f} + \textbf{g}_{\mathtt{P}} u_{\mathtt{P}}^*, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t))$ are linearly independent. Thus, the control inputs $u_{\mathtt{V}}^*(t)$ and $u_{\mathtt{P}}^*(t)$ can not be simultaneously singular on $I$. This concludes our proof. □

4.2.   Non-singularity of the optimal control input $u_{\mathtt{V}}^*$

First, we redefine (4.4) in terms of $u_{\mathtt{V}}^*(t)$. By using the relation $[\textbf{f} + \textbf{g}_{\mathtt{P}} u_{\mathtt{P}}^* + \textbf{g}_{\mathtt{V}} u_{\mathtt{V}}^*, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t)) = [ \textbf{g}_{\mathtt{P}}, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t)) u_{\mathtt{P}}^*(t)$ (since $[\textbf{f}, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t)) = 0$ and $[ \textbf{g}_{\mathtt{V}}, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t)) = 0$), we obtain the following:

where

Before we state our main result, we state the following assumptions that will be essential for what follows. It is important to note that the weights $c_{\mathtt{P}}, c_{\mathtt{V}}$, and thresholds $u_{\mathtt{Pmin}}, u_{\mathtt{Vmax}}$ are parameters that policy makers are free to decide at the onset of a pandemic. We enforce the effective cost of protection $c_{\mathtt{P}} (1 - u_{\mathtt{Pmin}})$ to be lower than the infection cost $c_{\mathtt{I}}$ to incentivize the protection adoption. In addition, recall that we restrict our analysis to the class of immunocompromised diseases such that $\hat{\beta} > \beta$ holds. These assumptions motivate the following mathematical conditions in the form of Assumptions 4.2. A further discussion on the choice of parameters is included in Section 5.

Assumption 4.2. We assume that the weighing and model parameters satisfy the following inequalities:

(i) $(\beta - \hat{\beta}) (c_{\mathtt{I}} - c_{\mathtt{P}}) \neq -\beta c_{\mathtt{V}} \gamma$;

(ii) $(\beta - \hat{\beta}) c_{\mathtt{I}} + \beta \hat{\beta} c_{\mathtt{V}} - \beta c_{\mathtt{V}} \gamma < 0$;

(iii) $(\beta - \hat{\beta}) (c_{\mathtt{I}} - c_{\mathtt{P}} (1 - u_{\mathtt{{Pmin}}})) + \beta \hat{\beta} c_{\mathtt{V}} u_{\mathtt{{Pmin}}} < 0$.

Assumptions 4.1 and 4.2 are sufficient conditions under which Theorem 4.1 holds.

Theorem 4.1. Suppose Assumptions 4.1 and 4.2 hold; then, the candidate optimal control input $u_{\mathtt{V}}^*$ cannot be singular on any open time-interval $I \subset [0, T]$.

Proof. The proof is by contradiction. Assume that the control input $u_{\mathtt{V}}^*$ is singular on $I$. A singularity in $u_{\mathtt{V}}^*$ is obtained when the the switching function $\phi_{\mathtt{V}}$ vanishes over $I$, which, in turn, implies $\phi_{\mathtt{V}}(t) = \dot{\phi}_{\mathtt{V}}(t) = \ddot{\phi}_{\mathtt{V}}(t) = 0$ for every/almost every $t \in I$. Equating $\dot{\phi}_{\mathtt{V}}(t)$ in (4.5) to zero and leveraging the fact that $[\textbf{f}, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t)) = \textbf{0}$, implies the following:

As a result of (4.10), the second time-derivative of the switching function in (4.9) is given by the following:

As the control input $u_{\mathtt{P}}^*(t) \neq 0$ (since $u_{\mathtt{P}}^*(t)$ lies in the interval $u_{\mathtt{P}}^*(t) \in [u_{\texttt{{Pmin}}}, 1]$, with $u_{\texttt{{Pmin}}} > 0$), on equating $\ddot{\phi}_{\mathtt{V}}(t) = 0$, we obtain the following:

Now, we express the vector $[ \textbf{g}_{\mathtt{P}}, [ \textbf{g}_{\mathtt{P}}, \textbf{g}_{\mathtt{V}}]](\textbf{z}^*(t))$ as follows:

where $\epsilon$, $\mu$, $\kappa: \mathbb{R}^3 \rightarrow \mathbb{R}$ are given by (dependency on time has been dropped for clarity)

It can be shown that the above obtained functions $\epsilon$, $\mu$, and $\kappa$ are indeed unique, since the three vectors $\textbf{g}_{\mathtt{V}}(\textbf{z}^*(t))$, $[ \textbf{g}_{\mathtt{P}}, \textbf{g}_{\mathtt{V}}](\textbf{z}^*(t))$, and $[\textbf{f}, [ \textbf{g}_{\mathtt{P}}, \textbf{g}_{\mathtt{V}}]](\textbf{z}^*(t))$ form a basis of $\mathbb{R}^{3}$ for each $t \in I$. Note that under Assumption 4.2(i), when $(\beta - \hat{\beta}) (c_{\mathtt{I}} - c_{\mathtt{P}}) \neq -\beta c_{\mathtt{V}} \gamma$ holds, $\kappa$ is well-defined.

Rewriting the denominator of $u_{\mathtt{P}}^*(t)$ in (4.11) (i.e., $\langle \pmb\lambda^*(t), [ \textbf{g}_{\mathtt{P}}, [ \textbf{g}_{\mathtt{P}}, \textbf{g}_{\mathtt{V}}]](\textbf{z}^*(t)) \rangle$) in terms of the right-hand side of (4.12), we obtain the following:

Substituting (4.13) in (4.11) under the conditions $\phi_{\mathtt{V}}(t) = \dot{\phi}_{\mathtt{V}}(t) = \ddot{\phi}_{\mathtt{V}}(t) = 0$ results in $u_{\mathtt{P}}^*(t) = -\frac{1}{\kappa(\textbf{z}^*(t))}$. Since, we assumed $u_{\mathtt{V}}^*$ to be singular on $I$, as a consequence of Proposition 4.1, it follows that $u_{\mathtt{P}}^*$ must exhibit a non-singularity on $I$. In other words, for a singularity of the optimal control input $u_{\mathtt{V}}^*$ to exist, it is necessary that $u_{\mathtt{P}}^*$ is a bang-bang control (i.e., either $u_{\mathtt{P}}^*(t) = 1$, or $u_{\mathtt{P}}^*(t) = u_{\texttt{{Pmin}}}$ for every/almost every $t \in I$). When the control input $u_{\mathtt{P}}^*(t) = -\frac{1}{\kappa(\textbf{z}^*(t))} = 1$, the following holds:

Note that the right-hand side of (4.14) lies in the range of $[0, 2]$, whereas the left-hand side is strictly negative by Assumption 4.2(ii). Thus, (4.14) can not be true. Now, we consider the case when $u_{\mathtt{P}}^*(t) = -\frac{1}{\kappa(\textbf{z}^*(t))} = u_{\mathtt{{Pmin}}}$, where the following holds:

Again, the right-hand side of (4.15) lies in the range of $[0, 2]$, whereas the left-hand side is strictly negative by Assumption 4.2(iii). Thus, the structure of the singular control input $u_{\mathtt{P}}^*(t) = - \frac{\langle \pmb\lambda^*(t), [\textbf{f}, [ \textbf{g}_{\mathtt{P}}, \textbf{g}_{\mathtt{V}}]](\textbf{z}^*(t)) \rangle}{\langle \pmb\lambda^*(t), [ \textbf{g}_{\mathtt{P}}, [ \textbf{g}_{\mathtt{P}}, \textbf{g}_{\mathtt{V}}]](\textbf{z}^*(t)) \rangle}$ does not hold, which implies that $\ddot{\phi}_{\mathtt{V}}(t) = 0$ is not possible on $I$. By integrating $\ddot{\phi}_{\mathtt{V}}(t)$ twice, it is deduced that ${\phi}_{\mathtt{V}}(t) = 0$ on $I$ is also not possible. Thus, under Assumptions 4.1 and 4.2, the optimal control input $u_{\mathtt{V}}^*$ is non-singular. This concludes our proof on the non-singular behavior of $u_{\mathtt{V}}^*$. □

Note that the above result which characterizes the behavior of $u_{\mathtt{V}}^*$ as a non-singular input is based on Assumptions 4.1 and 4.2. The impact of relaxing these assumptions presents an interesting research avenue, which we plan to explore in the future.

4.3.   Practical implication of theory

Our theory demonstrated that there is no simultaneous singularity, nor is there any possibility of the vaccination input exhibiting a singularity under the stated Assumptions 4.1 and 4.2. This has important implications on the public health policy. Specifically, the non-singularity results guarantee that the optimal vaccination policy is the simplest possible bang-bang control, which is often considered a more practical and appropriate intervention in epidemiological settings (see e.g., ^[28,29]). The mathematical proof of Theorem 4.1 rules out the existence of singularities in the vaccination input, which align with this broader understanding.

Further work needs to be done to fully characterize the optimal bang-bang control policy (i.e., determination of treatment level and transition times between treatment and no-treatment). Implementing our proposed strategies in real-world scenarios may be challenging, as implementation often requires adherence by individuals to the prescribed policies, and it is often difficult to ensure full compliance from people. It is to be noted that the optimality of a non-singular control has only been proven for the class of diseases leading to a compromised immunity (i.e., the reinfection rate is higher than the initial infection rate). As part of future work, we plan to expand our analysis to include diseases that confer a partial immunity, for which $\hat{\beta}$ is less than $\beta$.

5.   Numerical results

Now, we illustrate the trajectories of the optimal control inputs $u_{\mathtt{V}}^*$ and $u_{\mathtt{P}}^*$, and the evolution of the SIRI dynamics through numerical simulations. We demonstrate different phenomena, through three different cases obtained by changing the parameters and costs. In the first two cases, we select the parameters and costs such that Assumptions 4.1 and 4.2 are satisfied. In the third case, we violate the assumptions and illustrate the presence of singularities. We choose $u_{\texttt{Pmin}} = 0.2$ and $u_{\texttt{Vmax}} = 0.9$. In addition, for the endemic equilibrium to exist, the chosen parameters also satisfy the inequality $\gamma < \hat{\beta} u_{\texttt{Pmin}}$. Accordingly, for the first two cases, we choose the following weighing and model parameters which satisfy the above mentioned assumptions, whereas for the third case, we choose $\hat{\beta} < \beta$, with $x_{\mathtt{S}}(0) = 0.8, \; x_{\mathtt{I}}(0) = 0.2, \text{ and } x_{\mathtt{R}}(0) = 0$, where $x_{\texttt{j}}(0)$ for $\texttt{j} \in \{\texttt{S}, \texttt{I}, \texttt{R}\}$ denotes the initial state for the susceptible, infected, and recovered fractions of the population, respectively. The different costs and disease parameters are included in Table 1. Note that a different set of parameters and costs will not violate the theoretical results proposed in Proposition 4.1 and Theorem 4.1, as long as Assumptions 4.1 and 4.2 hold. Thus, the main results remain robust to the choice of parameters, which are also illustrated in the numerical simulations. The choice of the reinfection and recovery rates are governed by the basic reproduction number, the commonly used as a metric in epidemiological studies, to determine the strength of an infectious disease. The disease spreads as long as the basic reproduction ratio is greater than one. Since we focus on the case of a compromised immunity, we choose an infection rate $\beta < \hat{\beta}$. Several studies emphasized the variability in the reproduction numbers of different COVID-19 viral strains. For example, the authors of ^[37] predicted the reproduction number to be around $2.2$ in the early phase of COVID-19 spread, whereas the authors in ^[38] estimated the Omicron reproduction number in certain parts of the world to be around $8$. Such a large variation is potentially due to the behavioral and environmental circumstances alongside viral mutations. We have chosen our infection and recovery rates which varies within a range of $R_{0} \in [1.32, 20]$. In the first case, we set the costs such that the protection cost, $c_{\mathtt{P}}$, dominates the other two costs, even though the effective cost of protection $c_{\mathtt{P}} (1 - u_{\texttt{Pmin}})$ is lower than the infection cost. In the second and third cases, the protection cost is the lowest. The infection cost $c_{\mathtt{I}}$ being highest ensures that individuals are incentivized to choose either protection or vaccination, thus reducing the infection spread. The initial conditions of the model reflect that, a large proportion of the population is susceptible to the disease at the start of the epidemic, while only a small fraction is initially infected. We set the initial conditions such that there are no recovered individuals at the beginning of the epidemic, which corresponds to a first wave epidemic.

Certain numerical methods exist to analyze the existence of singularities. Most of the existing methods are efficient on linear systems. Our system under study is non-linear, and such numerical methods require linearization of the system around the operating point. We numerically validate our analytical results using the numerical solver Quasi-Interpolation based Trajectory Optimization (QuITO), which is famous for solving constrained nonlinear optimal control problems (see ^[39]). QuITO uses a direct multiple shooting (DMS) technique to discretize the control trajectory into several segments, and then obtains the optimal solution by solving for the control inputs at the boundaries of these segments. The trajectories which correspond to the states and control inputs for the three cases of weighing and model parameters are illustrated in Figure 2.

The plots in the top row represent the control inputs, whereas the bottom panel illustrates the corresponding state trajectories. In the first two cases, when all assumptions are satisfied, we observe that the simultaneous singularity of control inputs, as well as the singularity in $u_{\mathtt{V}}^*$, are completely absent, thus validating Proposition 4.1 and Theorem 4.1.

In the first case (i.e., Figure 2a), when the protection cost is high, and when the infection prevalence becomes very low, we observe that the complete removal of protection is the optimal policy. This is illustrated by the switching in $u_{\mathtt{P}}^*$ from $u_{\mathtt{P}}^* = u_{\mathtt{Pmin}}$ to $u_{\mathtt{P}}^* = 1$ after 42 days. In the second case, when the protection cost is the cheapest, we observe a complete adoption of protection throughout the time-horizon, irrespective of the infection level. In the first two cases, an interesting observation is the behavior of the trajectory of the optimal control input $u_{\mathtt{V}}^*$. At first, the behavior seems counter-intuitive; even though a sufficient fraction of the population is susceptible, vaccination is not applied as an input. This is explained based on the infection and reinfection rates $\beta$ and $\hat{\beta}$, respectively. Recall that the parameters are such that the rate at which the susceptible agents get infected (i.e., $\beta$) is lower than the rate at which the recovered agents get reinfected (i.e., $\hat{\beta}$). The susceptible agents have two options available to them: either (incurring a cost $c_{\mathtt{V}}$) they transit to the recovered state ($\mathtt{R}$) by getting vaccinated where they are likely to get infected at a (comparatively higher) rate $\hat{\beta}$ or they remain in the susceptible state ($\mathtt{S}$) and get infected at a (comparatively lower) rate of $\beta$. Quite intuitively, the latter option seems optimal for the susceptible agents. In other words, the choice of applying vaccinations (and thereby incurring a vaccination weighing parameter), then transiting to the recovered state, and finally getting reinfected at a higher rate $\hat{\beta}$ is not optimal (note that the infection weighing parameter $c_{\mathtt{I}} = 5$ dominates in the current scenario). Instead, agents prefer to remain in the susceptible state by not getting vaccinated. Therefore, we observe the optimal vaccination control $u_{\mathtt{V}}^* \equiv 0$ being satisfied throughout the given time duration. This leads to the important observation that when the reinfection rate is higher than the initial infection rate in a SIRI model with a high infection cost, the susceptible individuals prefer to remain unvaccinated and get infected at a lower rate.

Simulations in Figure 2a and 2b confirm the absence of simultaneous singularities for both the inputs $u_{\mathtt{P}}^*$ and $u_{\mathtt{V}}^*$, and the non- singularity of $u_{\mathtt{V}}^*$. These findings confirm the theoretical results outlined in the paper. It is important to note that a constant input represents a special case of non-singular control. Our analytical results focus on the special case of compromised immunity, where the reinfection rate ($\hat{\beta}$) exceeds the initial infection rate ($\beta$). As previously discussed, due to the high reinfection rate, individuals find it optimal not to get vaccinated, thus resulting in $u_{\mathtt{V}}^* = 0$, which is the lowest possible limit.

Now we focus on the simulations obtained under the third set of parameters. We violate Assumptions 4.1 and 4.2 by selecting $c_{\mathtt{V}} = 3$, $\beta = 3$, and $\beta > \hat{\beta}$, which implies partial immunity. Figure 2c illustrates the smooth behavior of the optimal control, which implies a singular vaccination input $u_{\mathtt{V}}^*$, whereas the input $u_{\mathtt{P}}^*$ remains constant at its lower limit. Additionally, this implies that Assumptions 4.1 and 4.2 may be close to necessary for the existence of non-singular optimal control laws. A full analysis of the case in which $\beta > \hat{\beta}$ is an interesting area of research that remains to be explored in future work. The possible singularities in $u_{\mathtt{V}}^*$ motivate us to investigate the behavior of diseases with a partial immunity.

Before concluding this section, we summarize the main advantages of our analysis and techniques in a broader context. We have demonstrated the non-existence of singular control inputs by analyzing the values of the time-varying switching functions $\phi_{\mathtt{P}}(t)$ and $\phi_{\mathtt{V}}(t)$. Necessary conditions for a singularity require these switching functions, along with their higher derivatives, to vanish identically. Our methodology is robust and can be applied to any compartmental model where the dynamics include control inputs and the running cost is linear in these inputs. Although our results focus on the relatively less-explored SIRI reinfection model, the technique for determining whether the control inputs are singular or non-singular could be useful for other epidemic models with different forms of control inputs. Thus, our technique is not model-sensitive and remains effective across various systems. This wide applicability of our analysis is due to the fundamental concepts of vanishing switching functions and their higher derivatives, which do not depend on the specific details of the underlying model.

6.   Conclusions and future work

In this paper, we considered the problem of optimal vaccination and protection for the class of SIRI epidemiological models. The biological significance of our results lies in the establishment of sufficient conditions on the susceptibility to infection and reinfection, and the cost of prevention and vaccination. Specifically, the proposed SIRI model takes reinfection into account, which is a key characteristic of diseases that result in short-term immunity. During the COVID-19 pandemic, we observed that the reinfection rates, particularly due to variants such as Delta and Omicron, were higher than the initial infection rates. Similarly, other diseases also exist which impart a compromised immunity. The existence of such real-world infectious diseases justifies the focus of this work on a compromised immunity. We proved that it is impossible for both the optimal control inputs to be simultaneously singular, when the immunity is compromised. Then, we performed a detailed analysis on the existence of a singularity of the optimal vaccination control input, and obtained sufficient conditions under which singular arcs (for optimal vaccination control input) are suboptimal and a non-singular vaccination control is optimal. Additionally, it is important to note that bang-bang control is often considered a more appropriate intervention in practical epidemiological settings. The numerical results provided valuable insights into the optimal control structure and evolution of the epidemic under such control inputs. Additionally, we illustrated that higher reinfection rates render vaccinations ineffective as the control input.

Since our bang-bang control optimality results are only guaranteed in the regime of compromised immunity, an extension of the analysis for which the infection rate is higher than the reinfection rate (as also seen in Figure 2c) would be worthwhile. It will be worthwhile to extend our analysis to include diseases that impart a partial immunity to individuals. Establishing the existence of singularities, and deriving expressions of the singular controls remains as a future work. Furthermore, an empirical validation using real data (e.g., involving a controlled interventional challenge study) would be valuable. Another worthwhile extension of our work would be to include infection testing and contact tracing as additional control inputs when all states of the compartmental epidemiological model (e.g., the SAIRU model in ^[40]) are not directly observable.

Author contributions

Urmee Maitra: Conceptualization, formal analysis, investigation, methodology, validation, writing - original draft, review & editing; Ashish R. Hota: Conceptualization, formal analysis, writing - original draft, review & editing; Rohit Gupta: Conceptualization, formal analysis, writing - original draft, review & editing; Alfred O. Hero: Conceptualization, formal analysis, writing - original draft, review & editing. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

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

Acknowledgments

This research was partially supported by grant CCF-2246213 from the National Science Foundation. The authors thank anonymous reviewers for their helpful suggestions.

Conflict of interest

The authors declare no conflict of interest in this paper.