A control theory approach to optimal pandemic mitigation

The recent outbreak of the illness COVID-19, caused by the SARS-CoV-2 virus, has resulted in a pandemic with unprecedented impact on societies all over the globe. Mitigation measures included complete lockdowns of societal life, with severe psychic, social, and economic consequences et al. (2020); Enserink and Kupferschmidt (2020). Hence a major focus of governments is on designing containment strategies which are as mild as possible, but substantial enough to limit the severity of the outbreak in order not to overwhelm the health service system (HSS). This requires reliable forecast, based on careful collection of data on the fraction of infected citizens, as well as extensive simulation Enserink and Kupferschmidt (2020); Dehning et al. (2020).

We discuss the system in terms of a so-called SIR model Hethcote (2000), referring to the fraction of susceptible (S), infected (I), and recovered (R) citizens in the population. We identify the recovered with all those who are neither susceptible nor infected ($R=1-S-I$); the dynamics are thus fully described by a set of two equations:

$$\begin{array}[]{r c l}\partial_{t}S&=&-\beta SI\,,\\ &&\\ \partial_{t}I&=&\beta SI-\frac{I}{\tau}\,,\\ \end{array}$$

(1)

where $S,I\in[0,1]$ are the fraction of susceptible and infected individuals in the population, respectively. Note that $I(t)$ denotes the actual fraction of acutely infected citizens at time $t$, no matter whether or not the infection has been realized by the individual, or has even been recorded. The infection of a susceptible individual by an infected is described by the infection rate $\beta>0$, while $\tau$ is the average duration of the infection of an individual until her recovery.

The task we address in this study is to limit, during the whole period of the pandemic, the current number of infected individuals such as to prevent the number of those needing intensive care from exceeding the capacity of the deployed HSS. Such control may be described by a control parameter $\alpha(t)$, which quantifies the effect of mitigation strategies upon the infection rate. We may write $\beta=\beta_{0}(1-\alpha)$, such that $\alpha=0$ and $\alpha=1$ correspond to usual societal life and complete mutual isolation of citizens, respectively.

In order to define $\alpha$ in a general way, we state that a certain value of $\alpha=1-\beta/\beta_{0}$ denotes the subset of all possible mitigation measures which lead to an infection rate $\beta\leq\beta_{0}$. We thus do not need to refer to any specific measures, but can formulate our approach in a very general way. The (more or less) accurate determination of these subsets is then the task of careful social (e.g., infection history) data analysis among citizens. This is illustrated in Fig 1, where mitigation strategies, followed by the public authorities, are indicated by the dashed and dotted curves, within a space spanned by the effect of the measures upon the infection rate, $\alpha$, and the cost incurred for economy and society as a whole, $f(\alpha)$.

As indicated by the explicit time dependence of $\alpha(t)$, we follow a control theoretic approach. This is in some contrast to earlier treatments which have assumed mitigation measures to be constant over time Harko et al. (2014); Gros et al. (2020); Zlatic et al. (2020); Ferguson et al. (2020). Instead, we aim at determining the optimal function $\alpha(t)$ which minimizes the impact on society, while at the same time avoiding the HSS to become overloaded. At the end of the mitigation scenario, herd immunity shall be reached, so that the epidemic comes to an end without further control. We do not consider potential vaccination scenarios here (as is done elsewhere Perkins and Espana (2020); Djidjou-Demasse et al. (2020)), so immunization can only be achieved via infection with the virus in the present study.

A dimensionless quantity which is frequently used in epidemiology is the reproduction number, $R\coloneqq\beta\tau S$, which denotes the average number of susceptibles infected by one infected individual. At the beginning of an epidemic ($S\approx 1$) and without mitigation measures deployed ($\beta=\beta_{0}$), one observes the basic reproduction number $R_{0}=\beta_{0}\tau$. In case of the COVID-19 pandemic, typical estimates are $R_{0}\approx 3$ Liu et al. (2020) and $\tau\approx$ ten days et al. (2020); He et al. (2020). It has been shown before that the inherently random, network-like structure of interactions between individuals mainly results in a readjustment of $R_{0}$ Parshani et al. (2010). Hence we follow a mean field approach, disregarding small scale inhomogeneities of the system. We consider a homogeneous scenario, where $\beta(t)$ depends on time, but is spatially constant. Defining a normalized time variable, $\theta:=\beta t$, we can rewrite Eq (1) as

$$\begin{array}[]{r c l}\partial_{\theta}S&=&-SI\,,\\ &&\\ \partial_{\theta}I&=&SI-\frac{I}{R_{0}(1-\alpha)}\,.\\ \end{array}$$

(2)

The trajectory of the system in the $S$-$I$-plane, $I(S)$, can be obtained by dividing the equations displayed in Eq (2) by each other. This yields

$$\frac{dI}{dS}=\frac{1}{R_{0}(1-\alpha)S}-1\,,$$

(3)

which has the solution:

$$I(S)=\frac{\ln S}{R_{0}(1-\alpha)}-S+C\,,$$

(4)

where $C$ is a constant of integration. When no mitigation measures are in place ($\alpha=0$), we have $I(1)=0$ and thus obtain

$$I(S)=\frac{\ln S}{R_{0}}-S+1\,.$$

(5)

This is plotted as the dashed curve in Fig 2 for the case $R_{0}=3$. The maximum turns out to occur at $S_{\textrm{peak}}=1/R_{0}$, where it reaches a value of

$$I_{\textrm{peak}}=1-\frac{\ln R_{0}}{R_{0}}-\frac{1}{R_{0}}\,.$$

(6)

At $S=S_{\textrm{peak}}=1/R_{0}$, the population has reached herd immunity since from then on the number of infected citizens decreases until zero.

If the disease is serious, one is faced with the problem that with a fraction of $I_{\textrm{peak}}$ people being infected, the number of those in need of hospitalization or even intense care may exceed the capacity of the HSS. We denote by $I_{h}$ $<I_{\textrm{peak}}$ the maximum fraction of infected citizens which can be managed by the HSS. It is limited by infrastructural aspects, such as the availability of staff or the size and areal density of hospitals, and is indicated by the horizontal dotted line in Fig 2. Any substantial overshoot of the dashed curve over the dotted line constitutes a catastrophe, as a major fraction of the population will then not receive proper health care or treatment. This must clearly be avoided by means of suitable measures, such as reducing mutual contacts between individuals, banning major assemblies, reducing mobility etc., thus reducing the infection rate. Such measures are described by the parameter $\alpha(t)$, which is to be discussed next.

It is clear that the aforementioned measures will have a more or less substantial impact on society, mainly through their detrimental effects on economy, but also through other societal (e.g., cultural) damage. This may be described by means of a cost functional,

$$K\{\alpha\}\coloneqq\int f(\alpha(t))\,dt\,,$$

(7)

where the cost function $f$ corresponds to the mitigation strategy chosen, i.e., the curve chosen in Fig 1. It denotes the cost incurred at a given control $\alpha$, along with the assumption $\partial f/\partial\alpha\geq 0\ \forall\alpha$; later we will require $\frac{\partial^{2}f}{\partial\alpha^{2}}\geq 0$ ¹¹1Convexity of the cost function is a reasonable assumption since one can imagine that chosen mitigation measures become more and more costly as $\alpha$ approaches unity. ²²2In principle, $f(\cdot)$ could depend on time $t$ explicitly or involve memory, e.g., because prolonged measures cause disproportinally higher cost. For simplicity, we do not treat such cases here..

The control problem we choose to address is to find a control trajectory, denoted by the function $\alpha(t)$, such that the impact on society, as described by $K$, is being minimized under the constraint that $I(t)$ never exceeds $I_{h}$ (capacity of the HSS) and at the end of mitigation—at unknown terminal time $t_{e}$—herd immunity is reached, i.e., $S(t_{e})=1/R_{0}$ (end of phase II in Fig 2).

We thus need to

$$\begin{split}\textrm{minimize}\quad&K\{\alpha\}=\int_{0}^{t_{e}}f(\alpha(t))\,dt\,,\\ \textrm{such that}\quad&\dot{\bm{x}}(t)=\bm{h}(\bm{x},\alpha(t))\,,\,\bm{x}(0)=\bm{x}_{0}\,,\\ &I(t)\leq I_{h}\,,\,S(t_{e})=1/R_{0}\,,\end{split}$$

(8)

where $\bm{x}=(S,I)$, and $\bm{h}(\bm{x},\alpha(t))$ as given by Eq 1 (with $\beta=\beta_{0}(1-\alpha)$). Minimization of mitigation time is covered by setting $f(\cdot)\equiv\mathrm{const}$.

Solving Eq 8 can be recast into minimization of the following functional:

$$J\coloneqq\int_{0}^{t_{e}}f(\alpha(t))+\bm{\lambda}(t)\cdot[\dot{\bm{x}}(t)-\bm{h}(\bm{x},\alpha(t))]+\mu(t)(I-I_{h})\,dt\,,$$

(9)

where $\bm{\lambda}(t)=(\lambda_{S}(t),\lambda_{I}(t))$, $\mu(t))$ are Lagrange multipliers. The introduction of $\mu(t)$ for the inequality constraint introduces additional constraints on $\mu(t)$, namely $\mu(t)\geq 0$ and the complementary slackness condition $\mu(t)(I^{*}-I_{h})=0$. These are also known as KKT conditions Kuhn and Tucker (1951). The star (^∗) represents the optimal quantities. Additionally, $S(t_{e})=1/R_{0}$ and $\bm{x}(0)=\bm{x}_{0}$ need to be enforced.

The necessary conditions for optimality can be evaluated by setting the first variation of Eq 9 to zero (for a detailed derivation see appendix A), we obtain:

	$\displaystyle f(\alpha^{*}(t^{*}_{e}))-\bm{\lambda}(t^{*}_{e})\cdot\bm{h}(\bm{x}^{*}(t^{*}_{e}),\alpha^{*}(t^{*}_{e}))=0\,,$		(10)
	$\displaystyle\dot{\bm{\lambda}}(t)=-\bm{\lambda}(t)\cdot\bm{\nabla}_{\bm{x}}\bm{h}|_{\bm{x}^{*}(t)}+\mu(t)\,\bm{\nabla}_{\bm{x}}(I-I_{h})|_{\bm{x}^{*}(t)}\,,$		(11)
	$\displaystyle\partial_{\alpha}f|_{\alpha^{*}(t)}-\bm{\lambda}(t)\cdot\partial_{\alpha}\bm{h}|_{\alpha^{*}(t)}=0\,,$		(12)
	$\displaystyle\lambda_{I}(t^{*}_{e})=0\,,$		(13)
	$\displaystyle\mu(t)\geq 0\,,$		(14)
	$\displaystyle\mu(t)(I^{*}-I_{h})=0\,.$		(15)

In addition to these, one also has the optimal system dynamics $\dot{\bm{x}}^{*}(t)=\bm{h}(\bm{x}^{*},\alpha^{*}(t))$. The necessary conditions become sufficient conditions if $\bm{h}(\bm{x},\alpha)$ and $f(\alpha)$ are convex in $\bm{x}$ and $\alpha$ Mangasarian (1966). The former can be checked to be valid for the SIR model and the latter implies that $\frac{\partial^{2}f}{\partial\alpha^{2}}\geq 0$.

The task remains to find Lagrange multipliers ${\bm{\lambda}}(t)$ and $\mu(t)$ which simultaneously satisfy the above conditions. This task usually involves a numerical search for the initial conditions of the Lagrange multipliers and evolve the system of ODE’s until the terminal conditions given by Eqs 10 and 13 are met. We escape the numerical difficulties arising with this procedure by first guessing a solution and then finding the appropriate Lagrange multipliers which verify optimality.

If immune response acquired after recovery from an infection is permanent, the pandemic will last until herd immunity is reached at the end of phase II. This is when $S(t)=S_{1}=1/R_{0}$, as indicated by the left vertical dotted line in Fig 2. This is the start of phase III, in which the number of infected citizens decays with no mitigation measures anymore in place (i.e., at $\alpha=0$). We will now discuss the time we expect it to take until this point is reached. Using Eq (17) with $I_{c}=I_{h}$, we can write

$$dt=-\frac{\tau}{I_{h}}\,dS$$

(22)

and hence for the total duration of the pandemic, $T_{0}$, until herd immunity is reached,

$$T_{0}=-\frac{\tau}{I_{h}}\int\limits_{S_{0}}^{S_{1}}dS\approx\frac{\tau}{I_{h}}\left(1-\frac{1}{R_{0}}\right)\,.$$

(23)

Here we have exploited the fact that in all cases of practical relevance, $I_{h}$ will be very small as compared to unity. Consequently, the duration of phase I will be very small as compared to phase II, such that the evaluation of the true duration of phase I is of minor importance. As a very good approximation, we have simply set $S_{0}=1$ and neglected the impact of $\alpha(t)$ on the dynamics for the short period of phase I.

Let us finally consider a few scenarios for some typical parameters, as we have in the current COVID-19 pandemic. The maximum number of known acute infections in Germany in spring 2020 was around 100.000, which was well tolerable for the HSS. Depending upon the percentage of cases which are officially recorded, the actual number of infected citizens may be considerably larger. If we assume a factor of two here, corresponding to 200.000 cases, we have $I_{h}\approx 0.0025$ (given $\approx 80$ million citizens in Germany). On the other hand, if there are more, and if we also take into account that the HSS could well take a few more patients, we might also consider a scenario with 800.000 acute infections at a time, corresponding then to $I_{h}\approx 0.01$. In both cases we also have to vary the average lifetime of the immune state, $\rho$, and observe its effect on the duration of the pandemic.

The two sets of scenarios are represented in the graphs in Fig 6. The remaining fraction of susceptible citizens is shown as the solid curves, while the dashed curves represent the control parameter, $\alpha$. As before, we have assumed $R_{0}=3$ for the basic reproduction number. We take the value $\rho=931$ days mentioned above for another SARS-CoV strain as a reasonable estimate. Using $\tau=$ ten days, this corresponds to $\rho=93\tau$. In order to cover this case, we have used the values $\rho/\tau=\{50,93,200,\infty\}$. For Germany, an HSS capacity of $I_{h}=0.01$ (top (a) graph) would correspond to roughly $32\%$ of hospital beds ⁴⁴4We do not consider the need for intensive care units here, one reason being lack of knowledge about the fraction of ICU cases. ($500.000$ in total) utilized for patients with COVID-19, if we assume a hospitalization rate of $20\%$ ⁵⁵5This is a conservative estimate given that the number published by the Robert Koch Institute Robert Koch Institute (2020) (17%) is based on reported cases only; the true hospitalization rate might be significantly lower.. The bottom (b) graph corresponds to a smaller HSS capacity ($I_{h}=0.0025$), for Germany, corresponding to $8\%$ utilization of hospital beds.

From the steadily decreasing dashed curves representing $\alpha(t)$, it is obvious that the mitigation measures can be gradually alleviated as time proceeds. In the top (a) graph ($I_{h}=0.01$) for infinite immunity ($\rho\to\infty$), one would reach the end of mitigation measures after about two years ($=66.7\tau$, with $\tau=10$ days). This is, however, hardly realistic. For the more realistic case, $\rho/\tau=93$, it would take about three years ($\approx 114.9\tau$). For an HSS capacity of $I_{h}=0.0025$, bottom (b) graph, clearly, there would be no chance to ever reach herd immunity for $\rho/\tau=93$. Instead, one would not reach any further than $\alpha\approx 0.5$, which still corresponds to rather harsh measures.

It should finally be noted that the number of fatalities is limited for all cases where herd immunity is reached. In particular, if $\phi$ is the fraction of fatalities among those which are infected, the fraction of fatalities with respect to the population is $F=\phi(1-R_{0}^{-1})$ for infinite $\rho$. If $\rho$ is finite, we find

$$F=\frac{\phi T_{r}}{T_{0}}\left(1-\frac{1}{R_{0}}\right)=\left(1-\frac{1}{R_{0}}\right)\frac{\phi\ln(1-X)}{X}.$$

(31)

This is precisely the scaling described by the curve displayed in Fig 4a, which shows that already well below $X=1$, the death toll incurred by the herd immunity strategy becomes prohibitively high if immune response decay plays a significant role.

We have shown that for a wide class of cost functions, in order to reach herd immunity without vaccination, it provides an optimal control strategy to keep the effective reproduction number, $R$, at unity during the majority of the duration of the pandemic. Deviations which depend upon the specific form of the cost function are limited to a narrow time window and can be considered negligible for practical purposes.

Reducing $R$ can be achieved through various measures, e.g., increased hygiene, physical distancing, or contact tracing Bulchandani et al. (2020). Keeping $R$ at the critical value of unity—above which epidemic spreading sets in—is, however, hardly feasible in practice, due to uncertainties as well as observation delays concerning the effects of mitigation measures. Development of robust optimal control scenarios taking such uncertainties into account is left for future investigations.

In this study, costs incurred at time $t$ have been considered local in time. Cost functions nonlocal in time (with a memory kernel) would be an interesting extension but go beyond the scope of this work. Costs associated with the number of infections have not been considered explicitly. Instead we kept the number of infections below an upper bound, e.g., the capacity limit of the health service system (HSS). Of course there are societal costs due to infections even below the limit of the HSS. However, if herd immunity is the goal and vaccination is not available, then there is no way around a fraction of $1-1/R_{0}$ of the population going through the infection. Moreover, the effectiveness of specific mitigation measures can depend on the number of infections; while contact tracing is an efficient measure for low case numbers, local health authorities can be overwhelmed if case numbers are high Contreras et al. (2020, 2021). Therefore, the socio-economic costs for establishing a given reproduction number $R$ might depend on $I$ as well. Here we focused on simple costs functions as a starting point allowing for analytical treatment.

Explicit expressions for the expected duration of the pandemic have been given, and we have seen that the duration of the pandemic increases strongly as the average lifetime of the immune state decreases. In particular, we can conclude that in case the immune response to SARS-CoV2 decays in a similar manner as for the formerly encountered SARS-CoV1 strain Wu (2007), using infection mediated herd immunity as a vaccination strategy for SARS-CoV2 would require a substantial fraction of health system capacity dedicated to COVID-19 patients (see Fig 6). However, as a consequence of global mobility there may be more pandemics coming which show different infection and immune response behavior. We therefore think that our results should be borne in mind for future use, as they are of rather general nature.