Introduction

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 [1, 2]. 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 [2, 3].

We discuss the system in terms of a so-called SIR model [4], 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:

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 α(t), which quantifies the effect of mitigation strategies upon the infection rate. We may write β = β₀(1 − α), such that α = 0 and α = 1 correspond to usual societal life and complete mutual isolation of citizens, respectively.

In order to define α in a general way, we state that a certain value of α = 1 − β/β₀ denotes the subset of all possible mitigation measures which lead to an infection rate β ≤ β₀. 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, α, and the cost incurred for economy and society as a whole, f(α).

Fig 1

Space of mitigation measures.

Sketch of the space of possible mitigation measures (grey shade), spanned by their effect on the infection rate, α, and their socio-economic cost, f(α). Normal societal life is at the origin, while the upper right corner corresponds to total mutual isolation of all citizens, which is the strongest possible intervention. The dashed and dotted curves depict possible choices for mitigation measures. Such curves correspond to the cost functions referred to in the manuscript (see Eq (7).

As indicated by the explicit time dependence of α(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 [5–8]. Instead, we aim at determining the optimal function α(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 [9, 10]), 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 ≔ βτS, which denotes the average number of susceptibles infected by one infected individual. At the beginning of an epidemic (S ≈ 1) and without mitigation measures deployed (β = β₀), one observes the basic reproduction number R₀ = β₀τ. In case of the COVID-19 pandemic, typical estimates are R₀ ≈ 3 [11] and τ ≈ ten days [1, 12]. It has been shown before that the inherently random, network-like structure of interactions between individuals mainly results in a readjustment of R₀ [13]. Hence we follow a mean field approach, disregarding small scale inhomogeneities of the system. We consider a homogeneous scenario, where β(t) depends on time, but is spatially constant. Defining a normalized time variable, θ ≔ βt, we can rewrite Eq (1) as

Fig 2

Trajectories in the (S, I)-plane.

Dashed curve: trajectory with no mitigation starting at (S, I) = (1, 0), R₀ = 3. Horizontal dashed line: maximum load of the HSS, I_h (here we have set I_h = 0.1 for clarity, although this is unrealistically large). Solid curve: trajectory, (S*(t), I*(t)), for an optimal choice of α(t) (see Eq 16). The corresponding characteristic of α(S) follows the dash-dotted curve in phase II, the mitigation phase. There is no mitigation in phases I and III (α = 0).

If the disease is serious, one is faced with the problem that with a fraction of I_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_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 α(t), which is to be discussed next.

Optimal control problem

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,

The control problem we choose to address is to find a control trajectory, denoted by the function α(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₀ (end of phase II in Fig 2).

We thus need to

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

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

The task remains to find Lagrange multipliers λ(t) and μ(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.

Heuristic approach

Let us first consider what is necessary to keep the fraction of infected citizens at a constant value, I_c. Since S varies with time, dI/dt = 0 entails dI/dS = 0, and hence from Eq 3 we obtain

This is indicated by the dash-dotted curve in Fig 2. Note that α(t) does not depend on the value of I_c.

Next we consider the cost function for proceeding from some S = S₀ to some S = S₁ < S₀ while maintaining I = I_c. Inserting Eq 16 in Eq 1, we find

We use this substitution to express the cost function as

Hence if S₀ and S₁ are fixed, I_c must be as large as possible to minimize K. This now guides our heuristic: we should control α such as to maintain I = I_h for as long as possible.

If our guess is valid, the trajectory we have to follow in order to optimally control the pandemic is the one indicated as the solid curve in Fig 2. It starts at (S, I) = (1, 0) and proceeds until I = I_h is reached. This completes phase I of the process, during which we set α = 0. Mitigation measures are then deployed, such that α jumps upwards to the dash-dotted curve. It follows that curve all through phase II, hence keeping I = I_h constant. As S decreases, α is gradually reduced until it reaches zero at the end of phase II. All through phase III, α is maintained at zero, while I gradually decays to zero because R < 1. This ends the pandemic.

Validation of the solution

We now proceed to verify our heuristic solution. We focus on phase II, as this is where the pandemic will spend the most amount of time. To do this we ask the question: Is it true that if the pandemic starts with I₀ = I_h, then for all S₀ > 1/R₀ and for all the cost functions f(α(t)) such that $$ \frac{\partial f}{\partial \alpha },\frac{{\partial }^{2}f}{\partial {\alpha }^2}\ge 0$$, optimal pandemic control must keep I(t) = I_h until $$ S\left(t_e\right)=\frac{1}{R_0}$$ is reached?

As we will see, the answer to the above question is yes. We proceed by setting $$ {\alpha }^{*}\left(t\right)=1-\frac{1}{R_0{S}^{*}\left(t\right)}$$ and $$ S^{*}\left(t\right)=S_0-\frac{I_h}{\tau }t$$ for $$ t\in [0,t_e^{*}]$$ and $$ t_e^{*}=(S_0-\frac{1}{R_0})\frac{\tau }{I_0}$$. The terminal conditions for the dynamics are given by $$ {\mathit{x}}^{*}\left(t_e^{*}\right)=(\frac{1}{R_0},I_h)$$. We can substitute the terminal conditions in Eq 10 and get the terminal condition for λ_S (the terminal condition for λ_I is given by Eq 13). The task remains to find μ(t) such that Eqs 11 and 12 are satisfied simultaneously.

Let’s have a look at Eqs 11 and 12 after making the substitutions. We have

and

One can hence find the Lagrange parameter μ(t) as

Lastly, there is also the issue of non-negativity of μ. If we assume the convexity of f we have $$ \frac{{\partial }^{2}f}{\partial {\alpha }^2}\ge 0$$. Hence the second summand in Eq 21 is non-negative. $$ \frac{\partial f}{\partial \alpha }\ge 0$$ also implies that λ_S is monotonically increasing (Eqs 19 and 20). If the cost of zero control is zero, then the terminal condition Eq 10 implies that $$ {\lambda }_S\left(t_e^{*}\right)=0$$, thereby implying λ_S(t) ≤ 0 in the interval $$ [0,t_e^{*}]$$. This shows that for time independent cost functions f under the assumptions that $$ \frac{\partial f}{\partial \alpha },\frac{{\partial }^{2}f}{\partial {\alpha }^2}\ge 0$$ and f(0) = 0 our heuristic solution is optimal in phase II.

Numerical results

We have shown that an optimal trajectory starting on the boundary (I₀ = I_h) remains on that boundary. To obtain optimal control trajectories for arbitrary initial conditions, we perform direct numerical optimization using the software library PSOPT [16]. In Fig 3 we show the numerical solutions to the control problem Eq (8) for various cost functions f_i(α(t)) ∈ {α(t), α²(t), α³(t)}.

Fig 3

Numerical solutions for optimal control.

Optimal control trajectories for different cost functions f_i(α(t)) ∈ {α(t), α²(t), α³(t)}. The corresponding optimal terminal times, $$ t_e^{*}$$, are determined as {65.99τ, 66.13τ, 66.31τ}. I₀ = 0.0025, I_h = 0.01, R₀ = 3, $$ S\left(t_e^{*}\right)=R_0^{-1}$$. R_∞ is the asymptotic reproduction number for t → ∞, given by R_∞ = −W(exp(−1 − R₀ I_h)), with the Lambert W function.

Clearly, in all scenarios the optimal trajectory I* reaches the threshold value I_h and remains there until $$ S^{*}\left(t_e^{*}\right)$$ is reached (phase II). Phase I, however, depends on the cost function applied. For linear costs, α(t) = 0 until I = I_h. With higher order cost terms, we observe non-zero control from the very beginning (see Fig 3). This is to reduce the amount of time spent at large control values α and thereby the total integrated costs. The optimal terminal time $$ t_e^{*}$$ increases with the order of the cost function (see Fig 3). We should note, however, that the influence of the functional form of f(α(t)), as expressed in the different shapes of the numerically derived curves, is minute, since the time axis is logarithmic, and the deviations are noticeable only during a very small fraction of time. Hence we see that the influence of the cost function, which corresponds to the chosen mitigation strategy, is finite, but can be regarded as negligible for practical purposes.

Duration of the pandemic

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/R₀, 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 α = 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

Influence of immune response decay

The introductory discussion was based on the idea that recovered patients stay immune for all times. However, it is well known that for some diseases, in particular of the SARS-CoV type, the immune response tends to decay after some time [17]. Hence there is some finite probability that recovered patients become susceptible again.

We now assume that the transition from the recovered to the susceptible state can be described as a Poisson process. In other words, we assume the probability that a randomly chosen, formerly infected citizen becomes susceptible in a time interval [t, t + dt] to be proportional to dt and independent of t. This modifies the dynamical system (1) to

with β = β₀(1 − α), and ρ the average life time of the immune state, averaged over all formerly infected individuals. Note that we conceptually include those who fell victim to the disease and thus do not become susceptible again. Their contribution to the average resusceptibilization frequency is zero, which merely increases the average immune lifetime, ρ. From the data in [17], we find that after three years the average IgG immune response against SARS-CoV had decayed to 55.6 percent. For a corresponding Poissonian process we can estimate ρ ≈ 931 days.

In Eq 24, we see immediately that the conditions to fulfill ∂_tI = 0 have not changed with respect to Eq 1. Hence the optimal control trajectory still obeys $$ \alpha \left(t\right)=1-\frac{1}{R_{0}S\left(t\right)}$$. In phase II, with optimal control, we obtain

with the solution

Again, herd immunity (and hence the end of the pandemic) is reached when S = 1/R₀, at a time we call T_r, referring to resusceptibilization (i.e., decaying immune response). Inserting this into Eq 26 yields

and hence

Although this looks rather awkward, it can be rewritten conveniently in terms of the pandemic duration, T₀, which we would find for infinite ρ. Defining the variable X = T₀/(τ + ρ), we find

X is the total duration of the pandemic if no loss of immunity occurs, divided by the average time it takes for a patient from her infection to the loss of immunity after recovery, τ + ρ. If immunity lasts very much longer than T₀, X is small. In this case, the logarithm in Eq 29 can be expanded and we recover T_r ≈ T₀. If, however, ρ is of order T₀ (remember that T₀ ≫ τ in all relevant cases), T_r diverges. This behavior is summarized in Fig 4a, in which X is chosen as the abscissa. We see that the duration of the pandemic becomes uncomfortably large when the total time from infection to resusceptibilization, τ + ρ, comes close to the pandemic duration with infinite immunity, T₀ (vertical dotted line).

Fig 4

Duration of the pandemic and minimum health system capacity.

(a) The normalized duration of the pandemic, T_r/T₀, as a function of the variable X = T₀/(τ + ρ) (Eq 29). (b) Solid curves: The minimum required health system capacity $$ {\widehat{I}}_h$$ to reach herd immunity (Eq 30) as a function of the duration of immunity after recovery, for different values of R₀ (from 1.5 to 4.0 in steps of 0.5). Dotted curve: limit R₀ → ∞. Circles represent the scenario for ρ = 93τ. Open: I_h = 0.01. Closed: I_h = 0.0025.

We might now ask how many acute infections the health system must be able to deal with in order to reach herd immunity at all. This can be derived by demanding lim_t→∞ S(t) = 1/R₀ in Eq 26. It is readily shown that the health system capacity required for reaching herd immunity is given by

This is shown in Fig 4b for different values of R₀. The dotted curve represents the (unrealistic) limiting case R₀ → ∞.

Only with infinite immune response lifetime (ρ → ∞), we observe an exponential decay of I after herd immunity has been reached (see also Fig 3). To understand the long time dynamics after mitigation (phase III) for finite ρ, we draw the phase portrait (see Fig 5). There exist two fixed points, (I₁, S₁) = (0, 1), a saddle for R₀ ≥ 1, and $$ (I_{\infty },S_{\infty })=({\widehat{I}}_h,1/R_0)$$, a stable fixed point for R₀ ≥ 1 (see S1 Appendix for details on the linear stability analysis). So for any initial conditions with I₀ > 0, the uncontrolled system will approach the stationary state, (I_∞, S_∞). Interestingly, the stationary fraction of infected coincides with the minimal HSS capacity $$ {\widehat{I}}_h$$ needed to reach herd immunity.

Fig 5

Phase portrait of the uncontrolled SIR model.

Phase portrait, $$ (\dot{I}(S,I),\dot{S}(S,I))$$, for the uncontrolled SIR model (α = 0) with finite immune response (Eq 24). The solid green curve shows a trajectory in phase III, with initial conditions (circle) I₀ = I_h (capacity limit) and S₀ = 1/R₀ (herd immunity), The dashed curves (orange) show the nullclines, $$ \dot{I}=0$$ (for S = 1/R₀ or I = 0), and $$ \dot{S}=0$$ (for S = (1 − I)/(IR₀ ρ/τ + 1)). The stable fixed point (diamond) is given by $$ I_{\infty }={\widehat{I}}_h$$, S_∞ = 1/R₀. Parameters: R₀ = 3, ρ/τ = 93, I_h = 0.01.

We have thus shown that given $$ I_h>{\widehat{I}}_h$$, herd immunity can be reached in finite time during mitigation phase II. After mitigation measures have been released (phase III), I converges to $$ {\widehat{I}}_h$$ in the long time limit. Moreover, I remains below I_h (see Fig 5) since re-entering the regime I < I_h from above would require a trajectory to cross itself (not possible for an autonomous system of ODEs (Eq 24) with unique solutions).

A few scenarios

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 ≈ 0.0025 (given ≈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 ≈ 0.01. In both cases we also have to vary the average lifetime of the immune state, ρ, 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, α. As before, we have assumed R₀ = 3 for the basic reproduction number. We take the value ρ = 931 days mentioned above for another SARS-CoV strain as a reasonable estimate. Using τ = ten days, this corresponds to ρ = 93τ. In order to cover this case, we have used the values ρ/τ = {50, 93, 200, ∞}. For Germany, an HSS capacity of I_h = 0.01 (top (a) graph) would correspond to roughly 32% of hospital beds. We 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%. This is a conservative estimate given that the number published by the Robert Koch Institute [18] (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.

Fig 6

Typical pandemic scenarios for different average immunity loss times, ρ/τ ∈ {50, 93, 200, ∞}, corresponding to curves from right to left (or see color code), and different values for I_h, namely 0.01 in the top (a) graph and 0.0025 in the bottom (b) graph.

Solid curves: S(t). Dashed curves: α(t). The fraction of acutely infected citizens is kept at I_h in phase II until herd immunity is reached (S = 1/R₀, horizontal dashed line). If this is successful (if $$ I_h>{\widehat{I}}_h$$, see Eq 30) phase III begins: Mitigation measures are being released (α = 0) and S(t) oscillates around its limiting value S_∞ = 1/R₀. Other parameters: R₀ = 3, τ = 10 days.

From the steadily decreasing dashed curves representing α(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 (ρ → ∞), one would reach the end of mitigation measures after about two years (= 66.7τ, with τ = 10 days). This is, however, hardly realistic. For the more realistic case, ρ/τ = 93, it would take about three years (≈ 114.9τ). For an HSS capacity of I_h = 0.0025, bottom (b) graph, clearly, there would be no chance to ever reach herd immunity for ρ/τ = 93. Instead, one would not reach any further than α ≈ 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 ϕ is the fraction of fatalities among those which are infected, the fraction of fatalities with respect to the population is $$ F=\varphi (1-R_0^{-1})$$ for infinite ρ. If ρ is finite, we find

Conclusions

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 [19]. 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₀ 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 [20, 21]. 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 [17], 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.

Materials and methods

Simulation code and data can be found in the following repository: https://github.com/poss-group/covid19-control.