key: cord-0604219-oca9chag
authors: Chyba, Monique; Klotz, Taylor; Mileyko, Yuriy; Shanbrom, Corey
title: A look at endemic equilibria of compartmental epidemiological models and model control via vaccination and mitigation
date: 2021-09-03
journal: nan
DOI: nan
sha: d90ea8a5d0fb2af9791f7d469b1fa4d30b7822ef
doc_id: 604219
cord_uid: oca9chag

Compartmental models have long served as important tools in mathematical epidemiology, with their usefulness highlighted by the recent COVID-19 pandemic. However, most of the classical models fail to account for certain features of this disease and others like it, such as the ability of exposed individuals to recover without becoming infectious, or the possibility that asymptomatic individuals can indeed transmit the disease but at a lesser rate than the symptomatic. In the first part of this paper we propose two new compartmental epidemiological models and study their equilibria, obtaining an endemic threshold theorem for the first model. In the second part of the paper, we treat the second model as an affine control system with two controls: vaccination and mitigation. We show that this system is static feedback linearizable, present some simulations, and investigate of an optimal control version of the problem. We conclude with some open problems and ideas for future research.

When modeling epidemics, compartmental models are vital for studying infectious diseases by providing a way to analyze the dynamics of the disease spread over time. The most basic of such models is known as the SIR model, which groups the population into three compartments (susceptible, infectious, recovered) and has a simple flow where an individual moves from susceptible to infectious to recovered ( [5] ). An additional compartment called the exposed group, can be included to obtain the SEIR model ( [36] ). This model is better suited for diseases with a latent period, the time when an individual has contracted the disease but is unable to infect others.

Although the SEIR model is a more accurate portrayal of an infectious disease than the SIR model, as most infectious diseases have a latent period, one limitation of this model is that it assumes that when someone recovers from the disease, they are immune to it forever. This is unrealistic because one can lose immunity over time. The SEIRS model is used to rectify this issue as this model assumes that individuals in the recovered group are able to return to the susceptible group ( [6] ).

The COVID-19 pandemic highlighted a few shortcomings of the classical SEIR models, as the latter do not account for the possibility of individuals staying asymptomatic throughout the course of the infection ( [13, 29] ) but at the same time capable of infecting susceptible individuals ( [18, 37, 38] ). One of the first models that looked into such a possibility was due to E. Sontag and collaborators in [23] , where the authors studied the effects of social distancing. Here we generalize the standard SEIR model in a similar way by introducing the SE(R)IRS model (Section 2.1). As the acronym suggests, individuals in the E compartment can pass directly to the R compartment without ever entering the I compartment. We should point out that in this paper we think of the E compartment as representing infected but asymptomatic individuals, while the I compartment represents individuals who are both infected and symptomatic. We have chosen to retain the traditional labeling for simplicity, although some modern authors use the letter A to denote the asymptomatic compartment [23] . We also allow individuals from the E compartment to infect those who are susceptible, although at a lesser rate than those from the I compartment.

In Section 2.2 we generalize further, adding a V compartment representing vaccinated individuals. The resulting model, which we call SVE(R)IRS, is significantly more complicated. We are still able to characterize equilibria in terms of the basic reproductive number, but fail to prove an endemic threshold theorem; thus this section is an open invitation to future research.

Creating a model is one thing, while analyzing it is quite another. In the first half of the paper we have chosen to focus our analysis on the equilibria of the systems and their stability. This choice was also motivated by COVID-19 and what will be the "end" of the pandemic. As the concept of herd immunity has received much attention in both the media and academia ( [1, 2, 9, 16, 39] ), the same cannot be said about of the notion of endemic equilibrium. There is a mathematical foundation for the idea of herd immunity ( [27] ), but as we demonstrate in Section 4, this does not mean the disease is eradicated. It is compatible with what we consider the more relevant idea of an endemic equilibrium: that the disease will always exist (hopefully in small enough numbers to no longer characterize a pandemic). Moreover, the stability of such an endemic equilibrium would reflect the possibility that new outbreaks or variants could cause spikes in infections, but that over time these numbers would drift back towards some state of "new normal" (see [4] for a nice overview of stability for SIRS models). The main idea is to design maintenance strategies to control the dynamic of the spread of the disease (through a yearly vaccine or seasonal non-pharmaceutical measures) to stay in a neighborhood of a sustainable endemic equilibrium.

While this discussion makes clear that our models and objects of study are motivated by COVID-19, we hope that our contributions can be applied to other infectious diseases with similar characteristics, including those yet to be discovered.

As in any mathematical modeling, there is naturally a trade-off between a model's complexity and its accuracy. In many ways the classical SIR model is useful mainly due to its simplicity, making both mathematical analysis and simulations painless. But its accuracy may be consequently limited. On the other hand, much more complicated compartmental models, such as that in [15] , may represent the dynamics of the disease very well at the expense of being computationally difficult. We hope that, like the popular SEIRS model, the models introduced here strike a reasonable balance by being simple enough for elementary dynamical systems theory and computations, while proving more flexible and accurate than the SEIRS model. In particular, Theorem 2 characterizes the stability of both the endemic and disease-free equilibria in terms of the basic reproductive number R 0 using only basic theory. But ignoring the effects of vaccination greatly misrepresents the course of pandemics like COVID-19. Yet our SVE(R)IRS model, while more accurate, was just complicated enough that similar analyses failed and we could prove no such theorem. For these reasons we believe these models lie somewhere near the right balance of complexity and simplicity. To our knowledge, neither has appeared in the literature before.

The second half of the paper, Section 3, treats the SVE(R)IRS model as an affine control system with two controls: vaccination and mitigation. This approach is similar to that of E. Sontag and collaborators in such works as [24, 41, 40] . However, while the authors of the latter paper focus on a single input control system, with the control representing non-pharmaceutical mitigation measures, we consider a bi-input control system, with the second control representing a vaccination rate. We study the resulting control problem from several perspectives. In Section 3.1 we prove that the system is static feedback linearizable, giving the explicit feedback transformation, and observe that it maps equilibria to equilibria; we also analyze the boundary conditions of both the original and transformed systems. In Section 3.2 we choose conceptually realistic control curves and explore the corresponding trajectories in both the original and linearized systems. In Section 3.3 we fix one control, and consider the 1-input system from an optimal control perspective. For the time-minimal problem we employ the Maximum Principle and analyze the singular controls.

We should mention that applications of control systems in biological and medical fields has seen an immense contribution from E. Sontag. In addition to producing a plethora of very influential papers, he trained and mentored generations of talented scientists and mathematicians who are now continuing in his footsteps and work on further expanding the applicability of control systems in various scientific disciplines. 

• β is the transmission rate, the average rate at which an infected individual can infect a susceptible • n is the population size • 1/σ is the latency period • 1/γ is the symptomatic period • 1/ω is the period of immunity. All parameters are necessarily non-negative. Note that, for simplicity, here and throughout this paper we choose to present our models without vital dynamics (also known as demography), often represented by the natural birth and death rates Λ and µ. Also note that when ω = 0, this reduces to the usual SEIR model, and one can further recover the simple SIR model by letting σ → ∞.

As described in Section 1, we now think of individuals in the E compartment as infected but asymptomatic (some authors call this interpretation a SAIRS model), while individuals in the I compartment are infected and symptomatic. Therefore in our SE(R)IRS model certain asymptomatic individuals can recover without ever becoming symptomatic; the duration of the course of their infection is denoted by 1/δ. Moreover, asymptomatic individuals can indeed infect susceptible individuals, however they do so at a reduced rate when compared to symptomatic individuals; this reduction is accounted for by the parameter α. We may assume α ∈ [0, 1] and δ ≥ 0. Note that when α = δ = 0 we recover the SEIRS model.

The SE(R)IRS model can be visualized as

and the corresponding dynamical system is given by

Proposition 1. The SE(R)IRS basic reproductive number is

Proof. We follow [25] by computing R 0 as the spectral radius of the next generation matrix. First, it is easy to see that the system has a disease-free equilibrium at (S, E, I, R) = (n, 0, 0, 0). We compute F = αβ β 0 0 and V = σ + δ 0 −σ γ so our next generation matrix is

The basic reproductive number R 0 is the spectral radius of this operator, which is the largest eigenvalue αγ+σ δ+σ β γ .

If we use the fact that n = S + E + I + R is constant we can reduce (1)-(4) to a 3-dimensional system in S, E, I. The reduced equations are given by dS dt

In the sequel we will study this simpler version of the system, recovering the value of R when convenient.

Our first result, proved in the next two sections, is the following.

Theorem 2. If R 0 < 1 then the disease-free equilibrium is locally asymptotically stable and the endemic equilibrium is irrelevant. If R 0 > 1 then the endemic equilibrium is locally asymptotically stable and the disease-free equilibrium is unstable.

Proof. The theorem follows from Lemmas 3, and 4, and 5.

Here "irrelevant" means epidemiologically nonsensical, as certain compartments would contain negative numbers of people; it still exists mathematically. This theorem is sometimes known as the "endemic threshold property", where R 0 is considered a critical threshold. In [27] , Hethcote describes this property as "the usual behavior for an endemic model, in the sense that the disease dies out below the threshold, and the disease goes to a unique endemic equilibrium above the threshold." The SEIR version is derived nicely in Section 7.2 of [36] . The SEIRS version can be found in [35] .

2.1.1. Analysis of endemic equilibria. Our system (5)-(7) has a unique endemic equilibrium at

Note that this endemic equilibrium is only realistic if all coordinates are positive, which requires positive, which is equivalent to

In other words, we have the following.

Lemma 3. If R 0 > 1 then there exists a unique endemic equilibrium for the SE(R)IRS model. If R 0 < 1 then no endemic equilbrium exists.

The linearization of the reduced system at p is

As expected, this matrix does not depend on n. Unfortunately, the eigenvalues of L are not analytically computable for general parameters. Note that L is nonsingular for generic parameter values. But the determinant does indeed vanish if and only if R 0 = 1.

If R 0 > 1 then the endemic equilibrium is locally asymptotically stable.

Proof. We apply the criteria (12.21-12.23) from [17] to L. Note that the trace is obviously negative, so (12.22) is immediately satisfied. Now compute

which is clearly negative, so (12.21) is satisfied. Finally, we compute the bialternate sum of L with itself,

is also negative, satisfying (12.23). These values are somewhat realistic for COVID-19. While there are no known precise values, and many of these parameters depend on the specific virus variant, these are at least roughly in agreement with some of the literature. Specifically, we assume that an asymptomatic individual is 10% as infectious as a symptomatic one, that individuals who become symptomatic have seven day periods of latency and of symptoms, that individuals who never develop symptoms are infected for 14 days, and that the period of immunity is 90 days. We choose the population size of 100 simply so that compartment values can be interpreted as percentages of a generic population. When β = 0.4 we have The eigenvalues of L are

These three eigenvalues all have negative real part, so the equilibrium is stable. It appears to be a spiral sink, signifying epidemic waves ( [6] ), as shown in Figure 1 . 2.1.2. Analysis of disease-free equilibria. One easily checks that the system has a disease-free equilibrium at (S, E, I, R) = (n, 0, 0, 0). The linearization of the reduced system there is

This is singular if and only if R 0 = 1, just like L.

Lemma 5. The disease-free equilibrium is locally asymptotically stable if and only if R 0 < 1.

Proof. The eigenvalues of N are

It is not obvious, but algebra shows that all three eigenvalues are real since our parameters are positive: the discriminant simplifies to 4βσ + (αβ + γ − δ − σ) 2 . Now λ 1 is clearly always negative. Next, we have

Thus λ 2 is also always negative, and in fact we have

Finally, Mathematica shows that 

Here 1 − ρ represents the efficacy of the vaccine, 1/ψ is the duration of efficacy of the vaccine, and 1/φ is the rate at which people are vaccinated. The first two are intrinsic to the vaccine itself, while φ can be thought of as a control (see Section 3). We may assume ρ ∈ [0, 1] and φ, ψ > 0. The associated dynamics are given by:

The dynamics of this model are significantly more complicated than those of the SE(R)IRS model in the previous section. We prove the analogues of Lemmas 3 and 5. We were unable to prove the analogue of Lemma 4, although experimental evidence suggests that it does hold.

Proof. Again we follow [25] by computing R 0 as the spectral radius of the next generation matrix. We use tildes to not confuse with the compartment V . First, a short computation shows that the system has a disease-free equilibrium at p 1 = (S, E, I, R, V ) = ψn φ+ψ , 0, 0, 0, φn φ+ψ . We computẽ

The basic reproductive number is the spectral radius of this operator, which is the largest eigenvalue: β γ αγ+σ σ+δ ψ+ρφ ψ+φ . Here again, if we use the fact that n = S + E + I + R + V is constant we can reduce the system and work with a 4-dimensional system in S, E, I, V . The reduced equations are given by dS dt

In the sequel we will study this simpler version of the system, recovering the value of R when convenient.

2.2.1. Analysis of endemic equilibria. The SVE(R)IRS system has three equilibria. One of them, denoted p 1 , is the disease-free equilibrium that will be analyzed in Section 2.2.2. The other two, denoted p 2 , p 3 , correspond to endemic equilibria and are square-root conjugate to each other. However, p 2 has negative components and is thus biologically irrelevant, while p 3 has all components positive and thus represents a true endemic equilibrium. Formally, we have the following result, analogous to Lemma 3, whose proof appears in the Appendix.

Proposition 7. If R 0 > 1, then there exists a unique endemic equilibrium point for the SVE(R)IRS system.

The linearization of our system at p 3 is unwieldy; Mathematica cannot even determine when the determinant vanishes, let alone compute eigenvalues. It can, however, produce the characteristic polynomial, so the methods of [17] applied in the proof of Lemma 4 could potentially work. However, the characteristic polynomial is a complex expression and it is hard to tell whether the coefficients are positive; something similar is expected for G, the bialternate product of this matrix with itself. Therefore the stability of p 3 remains open. For now we limit ourselves to examples, which provide hope that the analogues of Lemma 4 may indeed hold for the SVE(R)IRS model.

(α, β, γ, δ, n, σ, ω, φ, ψ, ρ) = 1 10 , The four eigenvalues all have negative real parts, so the equilibrium is stable. Note, however, that they are not all real, in contrast to the disease-free equilibrium case (see Proposition 8 below). We again see the appearance of a spiral sink due to epidemic waves ( [6] ). See Figure 2 . Example 4. Again we keep the parameter values as in Example 2 but here we fix the initial condition (S, E, I, R, V ) = (90, 5, 5, 0, 0) and vary β over time. We begin with β = 0.5(R 0 ≈ 1.8) and run the dynamics until we approach the corresponding endemic equilibrium. We then bump β up to 0.75 and repeat, then finally bump β up to 1.785 and repeat again. This aligns with the observations in [15, 31] and the real world data displayed in Figure 4 . 

Analysis of disease-free equilibrium. Setting E = I = 0 and solving the resulting system yields the unique disease-free equilibrium

Linearizing at this point yields the matrix

As in the SE(R)IRS model, this matrix is singular if and only if R 0 = 1. Proof. The eigenvalues of N are

where large amounts of tedious algebra show that

This form of Z makes it apparent that Z is always positive and thus all four eigenvalues are always real. Now it is also clear that λ 1 and λ 2 are always negative. Next, we have

Thus λ 3 is always negative as well, and in fact we have

Finally, Mathematica shows that det N > 0 if and only if R 0 < 1. Since det N = λ 1 λ 2 λ 3 λ 4 and λ 1 , λ 2 , λ 3 < 0, this shows that λ 4 < 0 if and only if R 0 < 1. This proves the Proposition.

In this section we analyze the SVE(R)IRS model as an affine control system with two controls of the form:q (t) = F 0 (q(t)) + u 1 (t)F 1 (q(t)) + u 2 (t)F 2 (q(t)). (22) The vaccine can be thought of as a first control over the system to steer the variables to desired values. The second control represents non-pharmaceutical interventions such as lockdown, social distancing, and mask mandates, and can also account for virus mutations (and therefore a change in transmission rate) that impact the parameter β. Introducingq = (S, E, I, R, V ) t , we have:

where u 1 (t) = φ(t) represents the vaccination rate and u 2 (t) = β(t) represents the transmission rate.

3.1. Static feedback linearization. Given a control system, it may be possible to change the state and control variables in such a way that the new system is a linear control system. Historically, the first result in this direction is a theorem of Brunovský [12] which says that all controllable linear control systems may be put into a standard form, which is now famously known as Brunovský normal form. In differential geometry language, such normal forms are essentially contact distributions of mixed order, otherwise known as generalized Goursat bundles [42] [43] .

is called static feedback linearizable (SFL) if there is an invertible map (t, z, v) = (t, ϕ(x), ψ(x, u)) such that the given control system transforms to a Brunovský normal form dz dt = Az + Bv.

The first results concerning when a given nonlinear control system is SFL were given by Krener in [30] , as well as by Brockett in [10] , and then Jakubczyk and Respondek in [11] . Constructing explicit maps for SFL systems is harder, and that work was started by Hunt, Su, and Meyer in [28] , and then a more geometric approach based on symmetry was developed in [22] by Gardner, Shadwick, and Wilkens, and finally work of Gardner and Shadwick [21] , [19] , and [20] provided what is now known as the GS algorithm for static feedback linearization. The work of Vassiliou in [42] and [43] provides both a test for determining equivalence of a generalized Goursat bundles and gives a procedure to construct an appropriate diffeomorphism.

The bi-input control system (23) fails to be controllable; however, as mentioned before we may reduce by the total population requirement. Indeed, this yieldṡ

where q = (S, E, I, V ) t . Theorem 9. The control system (24) is controllable and SFL with Brunovský normal form given byż

and z = (z 0 , z 1 , z 2 , w 0 ) t with the new controls v 1 and v 2 .

Proof. It is easily checked that the control system (24) is bracket generating and hence controllable. While the authors implemented the procedure in [43] via MAPLE to determine the linearizing map Φ below, one could determine the map by careful inspection after noting that the non-drift part of equation (24) is linear in the state variables and control-affine in the control terms. The feedback transformation (z, v) = Φ(q, u) is given by

The inverse map is given by

+ nσ(δω − γδ − γσ)w 0 − n 2 σω(σ + γ + δ)),

where l = σ(σγ + δ(γ − ω)),

is the recovered population in terms of the z coordinates. Notice that the controls u 1 and u 2 as functions of z and v are written via the (S, E, I, V ) state variables, and u 1 (z) is written also in terms of u 2 (z). Explicit expressions in terms of z may be written, but are unnecessary.

The domain for the state variables q is [0, n] 4 ; however, we may shrink to a smaller box region Ω q to avoid the disease-free equilibrium. Since the map (25) is linear in the state variables z and must have differential of full rank, then the domain Ω q is mapped to a region Ω z which is a parallelepiped in the z coordinates.

We now wish to understand the possible values for the controls (v 1 , v 2 ). In the (S, E, I, V ) space, the controls (u 1 , u 2 ) take values in the rectangle U = [0, φ max ] × [β min , β max ], where 0 < β min < β max ≤ 1 and φ max ≤ 1. Indeed, this yields the following proposition which is helpful for analyzing optimal control problems, as it classifies the possible (v 1 , v 2 ) regions by vertex behavior. First, we introduce the following expressions for convenience:

where it is clear that N 2 (z), N 3 (z), D 1 (z) > 0 on Ω z and L(z, v 1 ) has positive slope on Ω z as a linear function of v 1 .

and Σ + Σ Σ − = Ω z where Σ + = z ∈ Ω z : φ max > ∆β (αE(z) + I(z)) n , Σ = z ∈ Ω z : φ max = ∆β (αE(z) + I(z)) n ,

with ∆β = β max − β min . Then U is mapped into one of three types of parallelograms in (v 1 , v 2 ) space: V + , V, or V − whose boundaries are

Proof. First, under the backward static feedback map (25), the u 1 (z, v) and u 2 (z, v) may be written in terms of the expressions in (26) as

The region U in combination with the above equations immediately yields the inequalities Figure 5 . Sample control curves u 1 (t) and u 2 (t).

In this hypothetical scenario, the epidemic begins with no vaccinations and a transmission rate of 1. At time 100, a lockdown measure forces the transmission rate down to 0.1. At time 200, vaccination begins and grows linearly, until time 300 when we have u 1 = 0.05. At this time of peak vaccination, the lockdown ends, and the transmission rate jumps back up to 0.8 (lower than at the start of the epidemic, but much higher than during the lockdown) and stays there until the final time 500. At time 300 the vaccination rate begins decreasing linearly until time 400, when it stabilizes at 0.025. See Figure 5 .

The corresponding value of R 0 begins very high, then drops below 1 during the lockdown and decreases further during the vaccination campaign, but creeps back above 1 as the lockdown ends and vaccination rates stabilize. See Figure 6 . In this figure we also display the five state variables over time for this particular control scenario with the initial condition (S, E, I, R, V ) = (90, 5, 5, 0, 0). Figure 6 . The values of R 0 and all five states (S, E, I, R, V ) plotted over time, for the controls in Figure 5 .

In Figure 7 we display the same trajectory but in the new coordinates from the feedback transformation Φ from Section 3.1. We show the states z 0 , z 1 , z 2 , w 0 and the controls v 1 , v 2 over time. Finally, note that in Figure 6 the state variables are trending toward the point p ≈ (17, 0.5, 0.5, 10, 71). This point is the endemic equilibrium corresponding to the fixed parameters and final values of the controls β = 0.8 and φ = 0.025. In general this is expected: if the controls tend towards constant values then the states will tend toward the endemic equilibrium corresponding to these constant values (as long as R 0 > 1). Thus in this scenario governments can control the shapes of the curves to potentially avoid overcrowding of hospitals or exhaustion of resources, but the long term trajectory of the epidemic depends only on the limiting values of the controls as t → ∞. In Table 1 we display some values of the endemic equilibrium for particular choices of these values. 3.3. Optimal control. In this section we initiate a discussion of the optimal control problem to bring the population to a desired state, such as decreasing the number of infections to an acceptable level as quickly as possible. We limit ourselves to show the non existence of singular trajectories. It is unclear how to prove it directly on the original bi-input system (24) , however using the fact that singular trajectories are feedback invariant and that they do not exist for the Brunovský normal form we can conclude. We consider an affine bi-input system of the forṁ q(t) = F 0 (q(t)) + u 1 (t)F 1 (q(t)) + u 2 (t)F 2 (q(t))

where q(t) ∈ R 4 .

3.3.1. Maximum Principle and Singular Extremals. Assume that there exists an admissible timeoptimal control u 1 : [0, T ] → [0, φ max ], such that the corresponding trajectory q(t) is a solution of equations (28) and steers the system from q(0) to q(T ). For the minimum time problem, the Maximum Principle, see [8] , implies that there exists an absolutely continuous vector p : [0, T ] → R 4 , p = 0 for all t, such that the following conditions hold almost everywhere:

where the Hamiltonian function H is given by H(q, p, u) = p, F 0 (q) + p, F 1 (q) u 1 . Furthermore, the maximum condition holds:

where U is the domain of control. It can be shown that H(q(.), p(.), u(.)) is constant along the solutions of (29) and is greater or equal to 0. A triple (q, p, u) which satisfies the Maximum Principle is called an extremal, and the vector function p(.) is called the adjoint vector.

For our system (24) , given that the domain of control is of the form U = [0, φ max ] × [β min , β max ], where 0 < β min < β max ≤ 1 and φ max ≤ 1, the maximization condition implies that time-optimal strategies are formed by a concatenation of bang and singular arcs. Following notations in [8] , singular trajectories satisfy for all t the constraint

This constraint defines a subset Σ 1 of the space R 4 × R 4 \{0}. In [7] , B. Bonnard discusses the relation between singular extremals and the feedback classification problem for nonlinear systems; one the main results shows that singular extremals are feedback invariant. Using Theorem 9 and the linearizing map Φ, we have shown that our system (24) is feedback equivalent to the Brunovský normal form. The latter is a linear controllable system and therefore has no singular extremals.

Proposition 12. For system (24) with the domain of control given by U as above, the timeoptimal control strategies are bang-bang.

Deriving the time-optimal synthesis is non-trivial due to the complexity of the image of the domain of control U under the map Φ; see Proposition 10. It could however be implemented numerically since the switching functions for the time-optimal problem associated to the Brunovský normal form with two inputs are simple. Moreover, it seems that other costs would be more interesting to study.

This paper touched upon two important aspects related to the evolution of a pandemic. First, the design and implementation of efficient strategies to control diseases spread in a population, which, as we saw during the COVID-19 pandemic, is of paramount importance. Second, understanding evolutionary scenarios toward a "new normal", which is key for decision and policy makers to take proper actions at the correct time and avoid a prolonged stay in the pandemic phase.

We discuss the above two aspects in more detail below, but we should also mention that this work introduced a large number of open questions as well as directions for generalization. The most obvious question is whether the endemic threshold property in Theorem 2 holds for the SVE(R)IRS model of Section 2.2. All the examples we have explored suggest that the answer to this questions is affirmative, but with so many free parameters an exhaustive search for counterexamples is challenging. A different set of open questions concerns generalizations of both our models and our results. In particular, a natural addition to either model would be the vital dynamics of birth and death rates, as is common in the SEIRS literature. Other generalizations could include nonlinear transmission or seasonal forcing, as well as the appearance of new variants which would lead to having coefficients such as β depend on time. An alternative research direction would be the study of global rather than local stability of equilibria. Global stability for classical compartmental models has an extensive literature; see the comprehensive review article [27] . In particular, global stability was treated in [33] and [35] for SEIR, and in [14] and [34] for SEIRS. The methods of these and related papers might also prove effective for our models.

Herd immunity versus endemic equilibrium. In addition to our perceived need to consider compartmental models more closely adapted to COVID-19 than the classical models, this work was motivated in part by our observation that the national dialogue concerning the post-pandemic future focused largely on the concept of herd immunity rather than that of endemic equilibria. According to [27] , one has herd immunity when the immune fraction of the population exceeds 1 − 1 R 0 . In this case, the disease "does not invade the population". However, for classical models as well as the models presented here, we know that as long as R 0 > 1 we still do not reach a disease-free equilibrium which would correspond to the complete eradication of the disease. To illustrate the shortcomings of the herd immunity concept, reconsider Example 3 which includes annual vaccinations. Then we have 1 − 1/R 0 ≈ 0.69. According to [27] , herd immunity thus occurs if more than 69% of the populations is immune to the disease. However, at the endemic equilibrium for these parameters, we have R ≈ 66 and V ≈ 7. Thus 73% of the population are immune, which is above the required threshold. But over 6% of the population are either in the E or I compartments (≈ 3.4% each). This shows that although we have technically achieved herd immunity, a large number of people still suffer from the disease. Ideally, sufficient vaccination could eradicate a disease completely. By Proposition 8, this could happen if we took φ large enough to force R 0 < 1, as the dynamics would trend toward the disease-free equilibrium. However, consider the following parameters, which are realistic for COVID-19: (α, β, γ, δ, n, σ, ω, ψ, ρ) = 1 10 , 3 10 , 1 7 , 1 14 , 100, 1 7 , 1 90 , 1 180 , 1 10 . Here we have left φ free as a control. We compute that R 0 < 1 if and only if φ > 1/282. Thus, in order to set a trajectory towards disease eradication, the population would need to be re-vaccinated more often than once per year, which seems unlikely as a long-term sustainable scenario. Therefore for these parameters it seems more realistic that the best we could hope for is an endemic equilibrium with relatively small portion of the population infected at any given time. With annual vaccinations we would have approximately 0.73% of the population in the E or I compartments.

Control Strategies. The control theoretic investigation initiated in Section 3 offers many potentially fruitful avenues for research. Proposition 10 provides a description of the corresponding domain of control in the Brunovský normal form, and can be used to analyze the switching functions for the static feedbcack linearization and therefore corresponding time-optimal controls. They can then be mapped back to the original systems via the inverse transformation. More importantly, there are a large number of alternative cost functions from which one could reformulate the optimal control problem of Section 3.3: in particular, one may want to minimize the number of infections, or amount of vaccines administered, or the time taken to reach an endemic equilibrium, all subject to various boundary and endpoint conditions. The idea is to translate these cost functions via the SFL transformation given explicitly in the proof of Proposition 9. The difficulty to solve these new optimal control problems expressed in the Brunovský normal form is to be determined, and it is expected that numerical techniques will come into play.

The third of the above equations readily yields I = (σ/γ)E. Plugging this into the rest of the equations and setting κ = σ+αγ nγ we obtain 0 = −βκSE + ω n − S − 1 + σ γ E − V − φS + ψV 0 = βκSE − (σ + δ)E + ρβκV E 0 = −ρβκV E + φS − ψV.

Note that E = 0 yields the disease-free equilibrium, so we may assume that E = 0. Adding the bottom two equations to the first one and dividing the second equation by βκE we get

From the first two equations we easily find that

Using the expressions for S and V in terms of E, we obtain the following quadratic equation: ρβκaE 2 + (ρφa + ψa − ρβκb)E − (φc + ψb) = 0, which yields the following two roots E +,− = 1 2ρβκa

− (ρφa + ψa − ρβκb)± (ρφa + ψa − ρβκb) 2 + 4ρβκa(φc + ψb) .

Note that φc + ψb = φ 1 − ρ ρn − n γ(σ + δ) β(σ + αγ) + ψ 1 − ρ n − n γ(σ + δ) β(σ + αγ) = = n 1 − ρ ρφ + ψ − (φ + ψ) γ(σ + δ) β(σ + αγ) = n(ρφ + ψ)

We thus see that if R 0 > 1 then both roots are real, and E + > 0 while E − < 0. The latter gives a biologically irrelevant state. Computing the values of S and V from E + we obtain

ρβκb + ρφa + ψa− (ρφa + ψa − ρβκb) 2 + 4ρβκa(φc + ψb) .

S + = −c + ρaE + = 1 2βκ

− 2βκc − (ρφa + ψa − ρβκb)+ (ρφa + ψa − ρβκb) 2 + 4ρβκa(φc + ψb) .

Rather than showing directly from these expressions that S + and V + are positive, we proceed as follows. Note that instead of expressing S and V in terms of E we could as easily express any two of these variables in terms of the remaining one. Hence, we can obtain quadratic equations for S and V . We don't need the whole equations, but we note that both of them have a positive factor in front of the quadratic term, while the free constant term in the equation for V is φ(ρb − c), and in the equation for S it is −(ρb − c)(βκc/a + ψ)/ρ. Now, ρb − c = σ+δ βκ > 0. Thus, the two roots of the equation for V have the same sign. Since E − < 0 when R 0 > 1, we see that one of such roots is b − aE − > 0, hence we also have V + = b − aE + > 0. Similarly, we note that if c ≥ 0 the two roots of the equation for S have opposite signs, and since one of them is −c + ρaE − < 0, we have S + = −c + ρaE + > 0. The latter inequality is obvious when c < 0.

Finally noting that I + = (σ/γ)E + > 0, we see that having R 0 > 1 yields a single biologically relevant endemic equilibrium.

Challenges in creating herd immunity to SARS-CoV-2 infection by mass vaccination

Five reasons why COVID herd immunity is probably impossible

Optimal control of a SIR epidemic with ICU constraints and target objectives

Stability analysis of an eight parameter SIR-type model including loss of immunity, and disease and vaccination fatalities

Modeling infectious epidemics

The SEIRS model for infectious disease dynamics

Feedback equivalence for nonlinear systems and the time optimal control problem

Singular Trajectories and their Role in Control Theory

A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2

Feedback invariants for nonlinear systems

On linearization of control systems, L'Académie Polonaise des

A classification of linear controllable systems

Estimating asymptomatic SARS-CoV-2 infections in a geographic area of low disease incidence

On the global stability of SEIRS models in epidemiology

COVID-19 heterogeneity in islands chain environment

COVID-19 herd immunity: where are we?

Conditions for a matrix to have only characteristic roots with negative real parts

Asymptomatic transmission, the Achilles' heel of current strategies to control Covid-19

Feedback equivalence for general control systems

An algorithm for feedback linearization

The GS algorithm for exact linearization to Brunovsky normal form

Feedback equivalence and symmetries of Brunowski normal forms

A novel COVID-19 epidemiological model with explicit susceptible and asymptomatic isolation compartments reveals unexpected consequences of timing social distancing

Minimizing the infected peak utilizing a single lockdown: a technical result regarding equal peaks, medRxiv

Perspectives on the basic reproductive ratio

First special section on systems and control research efforts against COVID-19 and future pandemics

The mathematics of infectious diseases

Design for multi-input nonlinear systems

SARS-CoV-2 transmission from people without COVID-19 symptoms

On the equivalence of control systems and linearization of nonlinear systems

A study of computational and conceptual complexities of compartment and agent based models, Networks and Heterogeneous Media

Mitigation planning and policies informed by COVID-19 modeling: A framework and case study of the state of Hawaii

Global stability for the SEIR model in epidemiology

Global stability of the SEIRS model in epidemiology

Dynamical behavior of epidemiological models with nonlinear incidence rates

What the data say about asymptomatic COVID infections

Asymptomatic transmission of covid-19

Herd Immunity: Understanding COVID-19

Universal features of epidemic models under social distancing guidelines

An explicit formula for minimizing the infected peak in an SIR epidemic model when using a fixed number of complete lockdowns

A constructive generalised Goursat normal form

Efficient construction of contact coordinates for partial prolongations

Acknowledgments. The first and third authors were supported by NSF grant 2030789. The first three authors gratefully acknowledge support from the Coronavirus State Fiscal Recovery Funds via the Governor's Office Hawaii Department of Defense.

which are precisely the intervals Int 1 (z) and Int 2 (z) respectively. Since N 2 (z), N 3 (z), and D 1 (z) are all positive on Ω z , we have that L(z, v 1 ) has positive slope as a linear function in v 1 and therefore the region Int 1 (z) × Int 2 (z) is a parallelogram for each z ∈ Ω z . Moreover, as z changes continuously in Ω z , each parallelogram is continuously deformed to a new parallelogram with same vertical edges and positively sloped v 2 bounds. The only major change in shape that is relevant for optimal control purposes is the position of the vertices of each parallelogram. Letbe the vertices of each parallelogram labeled counter-clockwise starting from the highest value of v 2 . Vertices η 1 z and η 3 z will remain the largest and smallest values of v 2 for each parallelogram for each z because of the positive slope condition on L(z, v 1 ); however, η 2 z and η 4 z are not fixed relative to each other as z varies. Thus, there are three possible parallelogram types,Thus, the case of η 4 z − η 2 z = 0 corresponds exactly to a hyperplane Σ ⊂ Ω z and Σ + and Σ − correspond to being on either side of this hyperplane such that η 4 z − η 2 z is positive or negative, respectively.It is worth explicitly stating that the inequalities defining Σ ± in Proposition 10 is a statement about the max vaccination rate compared to the total proportion of infectious individuals multiplied by the difference in transmissibility.Finally, computations in Mathematica provide the following desirable result.Proposition 11. The static feedback map Φ takes equilibria to equilibria. That is, the diseasefree equilibrium p 1 and the two endemic equilibria p 2 , p 3 in the S, E, I, R, V coordinates are all mapped to the equilibrium plane of the Goursat bundle, which consists of points of the formRemark. Alternatively, one can consider the single-input control system wherein the control is the vaccination rate only and the transmission rate is held fixed. In this case, the control system is also static feedback linearizable; however, despite the static feedback transformation being rational, the expressions are far more cumbersome. One can in principle recover this map from the bi-input system by fixing the control u 2 to be constant and then appending the resulting differential equation. The single-input SF map is available upon request.

Simulations. Here we provide an example of control curves and plot the corresponding trajectories in the original variables. We then transform both the states and controls according to the feedback transformation of Section 3.1. Here we take φ = u 1 and β = u 2 and all other parameters are as in Example 2.Appendix: Proof of Proposition 7To find endemic equilibria of the SVE(R)IRS system, we need to solve the following system: