A survey on Lyapunov functions for epidemic compartmental models

In this survey, we propose an overview on Lyapunov functions for a variety of compartmental models in epidemiology. We exhibit the most widely employed functions, and provide a commentary on their use. Our aim is to provide a comprehensive starting point to readers who are attempting to prove global stability of systems of ODEs. The focus is on mathematical epidemiology, however some of the functions and strategies presented in this paper can be adapted to a wider variety of models, such as prey–predator or rumor spreading.

An abundance of approaches and mathematical techniques have been employed to capture the many facets and details which describe the spread of an infectious disease in a population.
In particular, compartmental models remain one of the most widely employed approaches. In these models, a population is partitioned into compartments, characterizing each individual with respect to its current state in the epidemic. One can then write a system of Ordinary Differential Equations (from here onwards, ODEs) to study the evolution in time of the disease.
These models usually take their names from the compartments they consider, the most renowned one being the Susceptible-Infected-Recovered (SIR) model. The SIR models can be extended to SIRS models by considering the acquired immunity to be temporary rather than permanent, allowing Recovered individuals to become Susceptible again. Various compartments can be added, depending on the characteristic of the specific disease under study: Asymptomatic, Exposed, individuals going though a phase of Waning immunity and many others.
A remarkably useful tool for the study of this kind of models are Lyapunov functions, which ensure global (or, in some cases, local) asymptotic convergence towards one of the equilibria of the system.
Given a system of n ODEs X = f (X ) and an equilibrium point X * , we call a scalar function V ∈ C 1 (R n , R) a Lyapunov function if the following hold: 1. V attains its minimum at X = X * ; 2. V = ∇V · f < 0 for X = X * .
The classical definition of Lyapunov function requires also the conditions 3. X * = 0 and V (X * ) = 0; however, these amount to a change of coordinates in R n and a vertical translation of V , so we will accept the more general definition. The existence of such a function guarantees the global stability of the equilibrium X * , as orbits of the systems naturally evolve towards the minimum power level of V .
The Basic Reproduction Number R 0 is a well-known threshold in epidemics models. Usually, R 0 < 1 suggests Global Asymptotic Stability (from here onwards, GAS) of the Disease Free Equilibrium (from here onwards, DFE), whereas R 0 > 1 suggests GAS of the Endemic Equilibrium (from here onwards, EE). In more complex models, the aforementioned conditions on R 0 might not be sufficient to prove the GAS of either equilibria, especially in cases in which the EE is not unique. Lyapunov functions often explicitly involve R 0 to guarantee the extinction of the disease or its endemicity over time.
Unfortunately, given a generic system of ODEs, there is no universal way of deriving a Lyapunov function, nor to rule out the existence of one. Moreover, the higher the dimension of the system of ODEs is, the harder it usually is to construct a Lyapunov function. However, there exist a few Lyapunov functions which have proven quite effective in a variety of different models.
In this survey, we collect some of the most relevant functions available in the literature, to provide the reader with a series of options to apply to the model of their interest, depending on its formulation. We include an extensive bibliography to complement the essential information of each model we present. This will provide the reader with a convenient starting point to investigate the availability of a known Lyapunov function to analytically prove the asymptotic behaviour of their system of ODEs. For the sake of brevity, we do not repeat most of the proofs to show that the functions we present are, indeed, Lyapunov functions for the respective system of ODEs. The missing proofs can be found in the papers we cite when introducing each model. Consider a model with compartments X 1 , X 2 , . . . , X n . Then, the DFE has coordinates X i = 0 for all i ∈ I, where I is the set of the indexes of infectious compartments, and the EE, which we indicate with (X * 1 , X * 2 , . . . , X * n ), has all positive entries. It is important to remark that, in order to use the Lyapunov functions for the EE given below, one only need to know that the EE is unique; an exact formula for each entry X * i is not necessarily required. A vast majority of Lyapunov functions in epidemic modelling fall into one of the categories listed below.

Goh-Lotka-Volterra. The Lyapunov function for the EE when
for some constants c i ≥ 0 to be determined [6,8,10,11,13,[16][17][18][19][20][21][22][23][24][25][26]. These functions are adapted from a first integral of the notorious Lotka-Volterra prey-predator system, and were popularized by Bean-San Goh in a series of paper [27][28][29]. 3. Quadratic. The Lyapunov function for the EE when R 0 > 1 is of the common form for some constants c i ≥ 0 to be determined, or the composite form Some examples can be found in [17,[30][31][32][33]. 4. Integral Lyapunov. Lyapunov functions given as integrals over the dynamics of the model. The integration interval often start at some EE value X * i and ends at the same X i ; this construction is very convenient if uniqueness of the EE is guaranteed, but the exact values of the EE are hard (or impossible) to determine analytically. Integral Lyapunov functions are particularly useful when the model includes multiple stages of infection, and consequently the infectious period changes from an exponential distribution to a gamma distribution [12,[34][35][36][37][38]. Integral Lyapunov functions, albeit in different forms, are widely used in models, which incorporate explicit delay, such as systems of Delay Differential Equations (from here onwards, DDEs), and age-structured models. However, these fall beyond the scope of this paper, and we will briefly comment on them in Sect. 3. 5. Hybrid. A linear combination of the above, which often includes the Goh-Lotka-Volterra in at least a few of the compartments of the system [14,22,23,25,[39][40][41][42].
For some high-dimensional models, proving convergence to the EE might require additional tools, such as the geometric approach used in [2,23].
Lastly, we must notice that not all compartmental models only exhibit convergence to equilibrium. Some systems of autonomous ODEs may present stable or unstable limit cycles [43][44][45], homoclinic orbits [43] or even chaos [46]. In such cases, clearly, no global Lyapunov function may exist.
In the remainder of this survey, we will present various models and the corresponding Lyapunov functions, covering all the cases listed above.

Epidemic models
In this section, we present various compartmental epidemic models with the corresponding Lyapunov function(s). We present the models from the smallest to the largest, in terms of number of compartments. We refer to [47,48] for a basic introduction on compartmental epidemic models, and to [49] for a detailed exemplification of Lyapunov theory in this setting.
We provide a schematic representation of the flows in most of the systems we present. Flow diagrams can be useful to provide a visual, intuitive interpretation of the parameters involved in each system. Arrows between compartments indicate a change in the current state of individuals with respect to the ongoing epidemics, whereas arrows inward/outward the union of the compartments represent birth rate and death rate in the population. Often, these last two rates are considered to be equal, as this assumption allows the population to either remain constant or converge to a constant value, reducing the dimensionality of the system and (hopefully) its analytical complexity. However, some models include additional disease-induced mortality, to increase realism when modelling severe infectious diseases. We uniform the notation throughout the various models we present in this survey as much as possible, and provide a brief description of each parameter the first time it is encountered. We remark that each variable is assumed to be non-negative, since it represents a fraction of the population, but the biologically relevant region varies depending on the specific model we are describing.
Moreover, we illustrate the corresponding Lyapunov functions for 2D models, showcasing a selection of their power levels. The same procedure can be easily adapted to 3D models, but the corresponding visualizations can be hard to interpret in a static image.

SIS
The SIS model is characterized by the total absence of immunity after infection, i.e. the recovery from infection is followed by an instantaneous return to the susceptible class. The ODEs system, which describes this situation is where β is the transmission rate and γ is the recovery rate. Notice that the population N = S + I is constant, thus we can normalize it to N = 1. Moreover, since S + I = 1, we can reduce the system to one ODE, which involves only infectious individuals  [17], in which the authors consider a birth/immigration rate different from the natural death rate; moreover, they include an additional disease-induced death rate from infectious class. Thus, the population is not constant and the system of ODEs, which describe the model is where represents the birth/immigration rate, μ the natural death rate and δ the diseaseinduced mortality rate. System (2) always admits the DFE, In [17], a Lyapunov function for the DFE is defined as For simplicity, we show that V (S, I ) is a Lyapunov function only in the case = μ and δ = 0 (i.e., there is no additional disease mortality). In this case, assuming N (0) = 1, N (t) ≡ 1 for all t ≥ 0, and then differentiating (3) with respect to time, we obtain Instead, the Lyapunov function for the EE is built using a combination of the quadratic and logarithmic functions Again, considering = μ and δ = 0 and differentiating (4) with respect to time, we obtain From N = 0, we know that = μ(S * + I * ). Moreover, from I = 0 at the EE, we observe that β S * = μ + γ . Combining these two equalities, we can write The authors provide two more examples of Lyapunov functions for the EE, namely and Power levels of the functions (3), (4), (5) and (6) are visualized if Fig. 1. By definition of a Lyapunov functions, orbits of the corresponding system (2) evolve on decreasing power levels, and they tend to the corresponding equilibrium as t → +∞.
In [32] the author found a simpler Lyapunov function for the DFE when R 0 < 1, i.e.
However, this last Lyapunov function (7) only ensures that I → 0 as t → +∞. To complete the proof of the convergence of the system to the DFE, one needs in addiction to invoke LaSalle's theorem [50] (see also [51,Thm. 3.4]), as is indeed done in [32].
Considering the importance of this theorem, especially when combined with the use of Lyapunov functions, we include its statement here. In particular, this theorem implies that, if we can prove the approach of the disease to the manifold describing absence of infection and the uniqueness of the DFE, then the DFE is GAS.

SIR/SIRS
The SIR model is characterized by the total immunity after the infections, i.e. recovered individuals can not become susceptible again. A classical example for this scenario is measles. The ODEs system, which describes this situation is where β is the transmission rate and γ is the recovery rate.
If we assume that recovered individuals eventually lose their immunity, we obtain the SIRS model. Denoting by α the immunity loss rate, we obtain the following ODEs system It is clear that, if α = 0, system (9) coincides with system (8).
These models admit only the DFE; in order to have an EE, we need to add the demography to model (8) or (9).
In [17], the authors consider the following ODEs system System (10) admits the DFE, E 0 = (S 0 , 0, 0), and the EE, E * = (S * , I * , R * ), which exists if and only if R 0 = β μ(μ + γ + δ) > 1. In [17], the Lyapunov function for the DFE is defined as follows whereas the Lyapunov function for the EE is the combination of the composite quadratic, common quadratic and logarithmic functions as follows The authors also present other Lyapunov functions for SIR/SIRS models; in particular, they also cite [52,53], in which some variations of system (10) are showed. Other Lyapunov functions for SIR/SIRS epidemic models are in [49], in which the authors use a graphtheoretic approach.
In [32], the author proved that the quadratic Lyapunov function (7) of the SIS model applies to the SIR and the SIRS, as well.

SEIR/SEIS/SEIRS
In [6], the authors study both SEIR and SEIS models. Many real world examples present a phase of exposition to the disease, between susceptibility and infectiousness. The models presented thus far, albeit simpler to study, are unable to replicate this mechanism.
The authors first analyze a SEIR model with demography and constant population, in which the disease is transmitted both horizontally and vertically. Individuals infected vertically pass first in the exposed compartment. The ODEs system, which describe the model is and R = 1 − S − E − I . The vertical transmission of the disease is represented by the probabilities p and q of being born directly in the Exposed compartment, rather than in the Susceptible one, and is represented by the inward arrow in compartment E.
The authors first provide an equivalent system, performing the substitution (S, E, I ) −→ (P, E, I ), where P := S + p μ β . They then proceed to prove the GAS of the EE, using the following Lyapunov function Later, the authors analyze a situation, in which the recovery does not provide immunity, namely the SEIS model. They also assume that a fraction r of offspring of the infective hosts is born directly into the infective compartment. In this case, the ODEs system, which describes the model is and S + E + I = 1. Notice that, due to the population remaining constant in system (12), one could in principle reduce its dimensionality and consider it as a planar system. The authors prove the GAS of the EE using the following Lyapunov function A natural extension to these models is the SEIRS [2,54], in which one can combine the existence of an immune compartment and the loss of immunity. It is described by the following system of ODEs μR μI (13) where g ∈ C 3 (0, 1], g(0) = 0 (meaning, in absence of infectious individuals, the disease does not spread) and g(I ) > 0 for I > 0. The classical choice is g(I ) = I , as in systems (11) and (12). Assuming moreover the authors of [54] . They then prove GAS of the DFE of system (13) through the use of the following linear Lyapunov function whereas the GAS of the EE is proved with a more complex geometrical method in [2].

SAIR/SAIRS
One of the main challenges of the Covid-19 pandemic was the presence of asymptomatic individuals spreading the disease. Such individuals must clearly be somehow distinguished from symptomatic infectious individuals, as they are likely to behave like a susceptible individual. Even though their viral load, and hence infectiousness, might be smaller, they are more likely to get in close contact with susceptible individuals.
In [23], the authors consider a SAIRS model. The main difference between this kind of models and the SEIR is that both asymptomatic and symptomatic hosts may infect susceptible individuals. The immunity is not permanent, i.e. recovered individuals will become susceptible again after a certain period of time. Moreover, vaccination is included. The ODEs system, which describe this model is The global stability analysis of the EE has been performed for two variations of the original model, described in the following.
The first model analyzed is the SAIR model, i.e. the case, in which recovery from the disease grants permanent immunity. In this case, the corresponding Lyapunov function is the combination of the Lotka-Volterra Lyapunov functions for S, A and I where c 1 , c 2 > 0. The second model is the SAIRS model, with homogeneous disease transmission and recovery among A and I , i.e. β A = β I and δ A = δ I . In this case, it is possible to sum equations for A and I , defining M := A + I , reducing the dimensionality of the system. Thus, the Lyapunov function can be written as the combination of the square function and the Lotka-Volterra as follows where w > 0. The global stability in the most general case is proved similarly to [2].

More exotic compartmental models
The aforementioned models are some of the most commonly used in literature. In order to capture additional disease-specific nuances, these model can be modified or extended by adding new compartments. In this section, in particular, we present a more complex model, in order to showcase one example of the integral Lyapunov function. Some diseases, for example, present different stages of infection. In this case, an infected individual can progress between two or more stages before recovering. In [12], the authors perform the global stability analysis via an integral Lyapunov function of a general class of multistage models. In their model, infectious individual can move both forward and backward on the chain of stages, in order to incorporate both a natural disease progression and the amelioration due to the effects of treatments.
The system of ODEs, which describes the model is where θ(S) is the growth function, f (N ) n j=1 g j (S, I j ) is the incidence term, ζ i (I i ), 1 ≤ i ≤ n, denote the removal rates of the I i compartment. Moreover, for any i, j = 1, . . . , n, the functions φ i, j (I j ) represent the rate of the disease progression if i > j and the amelioration if i < j.
The corresponding Lyapunov function for the DFE is linear in the disease compartments, i.e.
where c 1 = R 0 and c i ≥ 0 for all i = 2, . . . , n. For the global stability of the EE the authors made some assumptions on the aforementioned functions. In particular, they consider the following integral Lyapunov function where τ, τ i > 0, for all i = 1, . . . , n. For a more in-depth explanation on the functions (·) and ψ i (·) we refer to [12,Sect. 5]. Diseases which present multiple virus strains, due to the existence of different serotypes of the virus or due to a mutation of the original disease, may need to be modelled differently. Dengue, tuberculosis and various sexually transmitted diseases are caused by more than one strain of a pathogen. Influenza type A viruses mutate constantly: an infection with one of its strains gives permanent immunity against that specific strain. However, the so called "antigenic drift" produces new virus strains, thus the hosts only acquire partial immunity, or no immunity at all. Modelling these types of diseases requires the inclusion of cross-protective effects, in which the immunity acquired towards one strain offers partial protection towards another strain based on their antigenic similarity. In [13], the authors consider an n strain model, both without immunity and with immunity for all the strains. Moreover, they analyze an MSIR model, in which the M compartment represents the proportion of newborns who possess temporary passive immunity due to protection from maternal antibodies. For all the three models, the authors use a linear Lyapunov function to prove the global stability of the DFE and a logarithmic Lyapunov function to prove the global stability of the EE.
Other compartmental models include e.g. control strategies. For new ongoing epidemics, the most immediate strategy is including quarantine and isolation of infectious individuals. For well-known epidemics, for which a vaccination is available, it is useful to incorporate a vaccinated individuals compartment V to keep track of the two possible immunities, disease and vaccine induced, respectively. Usually, vaccination does not confer permanent immunity, and after a certain disease-dependent period individuals become susceptible again. An example can be found in [14], in which the authors analyze a SIRV epidemic model with non-linear incidence rate. The global stability of the DFE is proved using as linear Lyapunov function the infectious compartment I and the global stability of the EE, instead, using a combination of a quadratic function in S and a logarithmic function in the compartments I and V .

Conclusion
In this survey, we presented the most widely used Lyapunov functions in the field of epidemic compartmental models. We focused on systems expressed as autonomous systems of ODEs. These models allow for various interesting generalizations, of which we provide a noncomprehensive list below.
One extension of the classic compartmental epidemic models is the so-called multi-group approach, see e.g. [35,55]. These models describe n communities, interacting with each other, and whose internal evolution follows a standard compartmental model. A first example of such a model is presented in [56], in which the authors consider a n groups SIS model. In order to prove the GAS of the EE, they use a results on Metzler matrices. In [49], the authors consider a heterogeneous SIS disease model, for which they provide Lyapunov functions both for the DFE and for the EE. For the latter, they use a complex graph-theoretic method, for the details of which we refer to the original paper. Global stability of EE via Lyapunov functions for multi-group generalization can be found also for the SIR [57], SIRS [58], SEIR [59] and SAIR/SAIRS model [24]. Notice that, due to the complexity of the models, some of them require additional technical assumptions to prove the global stability of the endemic equilibrium.
Other classes of models include interactions between human and vector population, i.e. animals, which transmit the disease to humans, or with the pathogens, such as viruses or bacteria. In both cases, authors often include a compartmental structure for the non-human population. Some examples of vector-host models are shown in [9,10,60]. Another example can be found in [30], in which a SIR-B compartmental model is considered. Here the "B" denotes the concentration of the pathogen in the environment.
All the models discussed thus far are described by only autonomous systems of ODEs. However, in order to increase realism, it is possible to use non-autonomous systems to describe the spread of an infectious disease. This is the case of systems, in which some parameters change in time [61,62], to describe seasonal changes, or in which the state variables depend on the previous state, i.e. the model includes a time delay [63,64]. In these cases, it is still possible to find Lyapunov functions to prove the global stability of the equilibria using other techniques, described for example in [50].
Another popular option is to explicitly include delay in the system, such as in [63,[65][66][67][68][69] and [25]. In the latter the authors perform the global stability analysis of a SEIQR model, in which Q denotes the quarantined individuals compartment. They explicitly include a latent period for the infection, transforming two of the ODEs in DDEs. The corresponding Lyapunov function includes the integration over an interval whose size is precisely the latent period.
Lastly, a widely adopted strategy is to explicitly include the "time since infection" [70][71][72][73][74] in age-structured models. This allows to explicitly take into account time heterogeneity in the spread of an infectious disease in a population. Other modelling techniques for epidemics not treated in this survey, which nevertheless allow for the existence of Lyapunov function, include fractional derivatives [75][76][77][78] and Stochastic Differential Equations (SDEs) [79][80][81][82]. These cases are outside of the scope of this project, and we leave them as inspiration for future works.
As a final remark, some recent results in a more theoretical approach to the topic are worth mentioning. They focus on existence and characteristics of such functions rather than on applications to epidemiological models. We refer the interested reader to [83][84][85][86][87][88] and the references therein.