key: cord-0455955-g2jn6qbh
authors: Ballesteros, Angel; Blasco, Alfonso; Gutierrez-Sagredo, Ivan
title: Exact closed-form solution of a modified SIR model
date: 2020-07-31
journal: nan
DOI: nan
sha: fba21fe428ec1ad69ee0054fb0a2a29995d07afa
doc_id: 455955
cord_uid: g2jn6qbh

The exact analytical solution in closed form of a modified SIR system is presented. This is, to the best of our knowledge, the first closed-form solution for a three-dimensional deterministic compartmental model of epidemics. In this dynamical system the populations S(t) and R(t) of susceptible and recovered individuals are found to be generalized logistic functions, while infective ones I(t) are given by a generalized logistic function times an exponential, all of them with the same characteristic time. The nonlinear dynamics of this modified SIR system is analyzed and the exact computation of some epidemiologically relevant quantities is performed. The main differences between this modified SIR model and original SIR one are presented and explained in terms of the zeroes of their respective conserved quantities. We recall that both models have been recently used in order to describe the essentials of the dynamics of the COVID-19 pandemic.

The dynamics of an epidemic outbreak is intrinsically nonlinear, and this nonlinearity is obviously reflected in any dynamical system amenable to model epidemics. The simplest among these models are the so-called deterministic compartmental ones, and a large number of them have been recently considered in relation to the COVID-19 pandemic (see, for example, [1, 2, 3] and references therein).

While exact solutions have been found for some two-dimensional compartmental systems, like for instance the SIS (susceptible-infective-susceptible) model [4, 5, 6] , this is by no means the case for three-dimensional ones. We stress that for an infectious disease that produces some kind of immunity, like COVID-19, two-dimensional models do not appropriately describe the dynamics since a third compartment is needed in order to account for the recovered individuals that do not become immediately susceptible.

Among three-dimensional models, the well-known SIR (susceptible-infective-recovered) model S = −β S I,İ = β S I − α I,Ṙ = αI,

proposed by Kermack and McKendrick [7] is probably the best known one. Despite its apparent simplicity, it has been succesfully used to predict relevant features of the dynamics of a number of epidemics, including the actual COVID-19 pandemic [8, 9] . Nevertheless, the system (1) does not admit an explicit exact analytic solution in closed form [10] (see [11, 12, 13, 14] for different approaches to the mathematical properties of the SIR solutions). Moreover, this lack of an explicit closed-form solution is not a unique feature of the SIR model, but seemed to be a common feature of three-dimensional (and higher dimensional) compartmental epidemiological models. Therefore, the study of their dynamics typically relies on techniques from dynamical systems theory and numerical studies. Although these techniques allow a deep understanding of their associated nonlinear dynamics, the simplicity and accurateness provided by exact simple solutions would be indeed helpful both from the mathematical and the epidemiological perspectives.

In this letter we present the exact analytical solution in closed form of the modified SIR systeṁ

where α, β ∈ R + . This system has been proposed [15, 16, 17] as a more realistic model than (1) when the recovered individuals are removed from the population (not only due to death, but also to quarantine or other reasons). The general solution of this modified SIR system is

This is, to the best of our knowledge, the first exact solution of a three-dimensional compartmental epidemiological model in closed form. We also analyze the model and show that, in the typical range of the model parameters, the dynamics of (2) is actually quite close to the one of the SIR model (1).

In the next Section we derive this exact solution by making use of the fact that any epidemiological threedimensional model has a conserved quantity, which in turn is straightforwardly derived from the the more general result (recently proved in [12] ) stating that any three-dimensional compartmental epidemiological model is a generalized Hamiltonian system. Moreover, the conserved quantity turns out to be just the Casimir of the Poisson algebra of the underlying Hamiltonian structure. In Section 3 we present the analysis of the modified SIR system (2) both from a dynamical systems approach and from a Poisson-algebraic point of view, and we show that the exact solution (3) is helpful in order to obtain some relevant epidemiological quantities in a simple and exact form. Finally, the main differences between the SIR and modified SIR systems are analysed in Section 4, where we show that these differences can be understood in terms of the the zeroes of their respective conserved quantities, which are again the Casimir functions for both models that are obtained through the formalism presented in [12] .

In order to find the exact solution of (2) we make use of the following recent result (see [12] for details). Proposition 1. [12] Every epidemiological compartmental model with constant population is a generalized Hamiltonian system, with Hamiltonian function H given by the total population.

For the system (2) the generalized Hamiltonian structure is thus explicitly provided by the Hamiltonian function

together with the associated Poisson structure, which is found to be given by the fundamental brackets

and leads to the system (2) through Hamilton's equationṡ 

Since every three-dimensional Poisson structure has a Casimir function C, i.e. a function C : U ⊆ R 3 → R such that {S, C} = {I, C} = {R, C} = 0, then C is a conserved quantity for any generalized Hamiltonian system (6) defined on such a Poisson manifold. Therefore: Corollary 1. Every three-dimensional epidemiological compartmental model with constant population has a conserved quantity, which is functionally independent of the Hamiltonian function.

Note that in case that H (4) is functionally dependent of C, the dynamics (6) would be trivial. For the specific Poisson algebra (5) the Casimir function is found to be

We can use this Casimir function to restrict the dynamics of (2) to the symplectic leaf defined by the value of C given by the initial conditions S(0) = S 0 , I(0) = I 0 , R(0) = R 0 , namely

This can be also used in order to reduce the system (2) to a nonlinear ODE, since from (7) and (8) we obtain the phase space equation

which can be inserted within (2) in order to get the following nonlinear ODE for the variable S:

This ODE suggests the change of variable

thus obtainingȦ

If we now set

The general solution to this ODE is a logistic function with characteristic time τ = (β − α) −1 , i.e.

The integration constant d is fixed by the initial condition B(0) = S0 S0+I0 , thus obtaining e d = I0 S0 . Therefore we can write

Now, inverting the change of variables (13) we get

and finally, from (11), we obtain

From the phase space equation (9) we directly get

and inserting (18) we are able to obtain I(t) without any further integration. Finally, we have that

Note that I(t) is related to S(t) by

and from the conservation of the total population, we find

Summarizing, equations (18) shows that the susceptible population follows a generalized logistic function, or Richards' curve, with characteristic time τ and the relevant constants set to satisfy that S(0) = S 0 and lim t→∞ S(t) = 0. Moreover, the dynamics of the infective population given by (20) is essentially this same function multiplied by an exponential with the same characteristic time. This is indeed a very natural dynamics for infective processes and, as we will see in the sequel, this dynamics strongly resembles the one described by the famous SIR model (1) , provided that the range of values for the parameters α and β is similar to the one found in actual epidemics.

Remark 1 It is worth stressing that the method here presented is indeed applicable to any three-dimensional compartmental model, provided we are able to find the Casimir function of the associated Poisson structure. Nevertheless, the distinctive feature of the system (2) is that the resulting ODE admits a closed-form solution. We recall that in [12] such Casimir function approach was used in order to find the solution for some epidemiological models in terms of an inverse function.

In this Section we briefly analyze the main features of the modified SIR system (2) . Without any loss of generality we can assume that R 0 = 0, so S 0 + I 0 = 1, and the solution of (2) reads

As we have previously stated, the behavior of S(t) is that of a generalized logistic function while the evolution of the infective population I(t) is given by a generalized logistic function times an exponential. This means that since βτ > 1, the logistic term dominates for large times and therefore lim t→∞ I(t) = 0. However during the first stage of the outbreak the exponential term is the dominating one, and thus the model presents the characteristic infection peak (for appropriate values of the parameters α and β). The behavior of the functions S(t) and I(t) for different values of α and β is shown in Figure 1 .

A fundamental question to be answered by any epidemiological model is whether, for given values of the parameters, there will be an outbreak. For the modified SIR system we see, simply by evaluating the second equation from (2) at t = 0, d dt t=0

and the outbreak will exist if and only if βS 0 > α(S 0 + I 0 ), or equivalently

which in the case S 0 + I 0 = 1 means that

Obviously, this same result can be obtained by checking the condition for which I(t) has a maximum. Moreover, the analytic solution allows us to exactly determine the time at which the infection peak t peak is reached, and we obtain

which is positive if and only if (25) holds. The fraction of infected population at the infection peak reads

Another interesting insight is gained by computing the area below the infective curve I(t). In order to do that, we do not even need to perform the integration of I(t), since from the third equation in (2) we get

This result is specially interesting from a parameter estimation point of view, since it allows to obtain a value for α directly from the data. Afterwards, assuming that S 0 and I 0 are known, β can be obtained, for instance, from (27). Thus, the exact solution in closed form greatly simplifies the fitting with actual data. Moreover, as we will see below, since the dynamics of the SIR and modified SIR systems are quite close (for a realistic range of the parameters), this procedure for the determination of the parameters of the modified model provides a good approximation for the parameters of the SIR one.

A related interesting quantity from the epidemiological point of view is the removal rate, defined byṘ(t) = αI. While for the SIR model it can only be approximated by a closed-form expression in certain limits (see [5] ), in the modified SIR system it can obviously computed exactly. Therefore, the behaviour of the removal rate (divided by α) for the modified SIR system can be directly extracted from Figure 1 .

For any epidemic outbreak, it is also enlightening to analyze the intersection of the susceptible S(t), infective I(t) and recovered R(t) functions. The closed-form solution of the modified SIR model allows us to get some exact results in this respect, which we write down in the following Proposition 2. For the modified SIR system given by (2), with β > α and initial conditions S(0) = S 0 , I(0) = I 0 , R(0) = 0 such that S 0 > I 0 > 0 and S 0 > α/β, any two of the curves S(t), I(t) and R(t) always intersect exactly once, regardless of the exact values of the initial conditions and parameters of the system. Furthermore:

i) The curves S(t) and I(t) intersect before the infection peak if β > 2α, exactly at the infection peak if β = 2α

and after the infection peak if β < 2α.

ii) The three curves S(t), I(t) and R(t) intersect in a common point if and only if β α < Proof. The solution of (2) when R(0) = 0 is given by (23). In particular, the unique time at which the curves S(t) and I(t) intersect can be explicitly computed from (21), and it reads

Since we are assuming S 0 > I 0 , this time is always positive.

From (23) we can also compute the times t SR and t IR such that R(t SR ) = S(t SR ) and R(t IR ) = I(t IR ). It is easy to check that these times are given by the common expression

where X is a solution of the equation

in the case of t SR , while X is a solution of the equation

in the case of t IR . An elementary computation shows that equation (32) has only one solution, and therefore t SR is unique. Equation (33) has 2 solutions, the first one living in (0, 1) and the second one within (1, ∞). The first of these solutions results in a negative or complex time, so t IR is defined by the unique solution of (33) in (1, ∞) . Therefore, we have proved that t SI , t SR and t IR are unique, which means that the curves S(t), I(t) and R(t) intersect exactly once.

To prove ii) we first note that t * = t SI if and only if X = 2S 0 , which is a solution of both (32) and (33) if and only if

Given that S 0 < 1, then 1 3 3 2 β/α < 1, and therefore β α < log 3 log(3/2) . Statement iii) is a consequence of i) and ii), since in order that the intersection coincides with the infection peak we need that β = 2α, and substituting this into the condition for triple intersection (34) we get S 0 = 3 4 . Equivalently, note that if β = 2α, then by (28), S(t peak ) = I(t peak ) = 1 4S0 , so R(t peak ) = 1 − 1 2S0 , and by imposing that they coincide we get S 0 = 3 4 . Remark 2 The condition S 0 > α/β in the previous Proposition is just the condition for the existence of an infectious peak (26) but we are not using it explicitly in the proof. Therefore, all the previous results which do not involve the infection peak are true if this condition is removed.

The fact that the functions S(t), I(t) and R(t) always intersect is a striking difference with respect to the original SIR system of Kermack and McKendrick (1).

Finally, it is worth performing a more detailed analysis of the dynamics of the modified SIR system (2) when compared to the original SIR model (1) . Although the exact closed-form solution here obtained in the modified case is valid for any values of α and β, for the sake of brevity from now on we only consider the case β > α (recall from (25) that this is the regime in which an actual outbreak does exist).

Qualitatively, the most significant difference between both dynamical systems is the stability of their fixed points. From equations (1) and (2) we see that I = 0 is a line of non-isolated fixed points for both models. The stability of any of these fixed points p = (S, 0) is defined simply by the trace of the Jacobian evaluated at this point, J(p). For the SIR system this trace is

while for the modified SIR system it reads

Therefore, for the modified SIR system all the points belonging to the line I = 0 (with the exception of (S, I) = (0, 0)) are unstable (recall that we only consider the epidemiologically relevant case β > α). Meanwhile, for the SIR system points such that S > α/β are unstable, while points with S < α/β are stable. This can be clearly appreciated in Figures 2 and 3 , where the corresponding flows are presented for both systems. Coloured curves correspond to the phase space equation I(S) for each model. Trajectories of the system starting at any point of the appropriate curve will follow this curve (in the direction of the flow) in order to reach the relevant fixed point.

The differences regarding the fixed point structure of these two systems can also be analyzed algebraically. For the SIR system (1), it is well-known that the phase space equation is

where C is a constant (this is the equation of the curves in Figure 2 , for different values of C). In fact,

is the Casimir function for the associated Poisson structure (see [12] for details). It is easy to prove that the equation I(S) = 0 always have a solution S ∈ (0, α/β). However, the phase space equation I(S) = 0 for the modified SIR system, where I(S) is given by (9) , always has S = 0 as a solution (see Figure 3 ). This equation is directly obtained from the Casimir function (7), and it is interesting to compare this Casimir function (7) with the exponential of (38). Remark 4 The previous discussion shows that the different qualitative behaviour of the systems (1) and (2) can be algebraically understood through the differences between the Casimir functions (7) and (38), and in particular, the different structure and location of their zeroes within the phase space.

From the epidemiological point of view, the existence of stable fixed points (different from (0, 0)) in the I = 0 axis explains the well-known fact that in the original SIR model of Kermack and McKendrick the whole population is not infected during the evolution of the infection. While these results can be obtained from a dynamical systems approach, it is interesting to note their direct connection with the algebraic and geometric structure of the Poisson manifold underlying their description of epidemiological models as generalized Hamiltonian systems. 

for any initial conditions such that I 0 = 0. This shows that the only stable fixed point is (S, I) = (0, 0). Note that for the modified SIR system it takes an infinite time to reach this fixed point, regardless of the initial conditions.

In Figure 1 , plots for some trajectories S(t) and I(t) contained in the phase space orbits from Figures 2 and 3 are depicted. In the left column β = 1 and α = 0.2 while in the right column β = 1 and α = 0.6. Each plot contains four different curves: coloured ones correspond to I(t) while black ones correspond to S(t), and solid ones correspond to the modified SIR system while dashed ones correspond to the original SIR system. The first row shows the dynamics for initial conditions such that the outbreak rapidly extinguishes. The second row shows the limiting case given by S 0 = α/β (note that this value is the same for both models since S 0 + I 0 = 1). The third row shows the typical behavior for values of the parameters and initial conditions for which there is an actual outbreak, and therefore I(t) has a maximum. The fourth row shows a situation such that at the beginning of the outbreak a fraction of the total population is immunized.

These plots show all possible qualitatively different dynamics for the SIR and modified SIR systems. Essentially, as far as the ratio α/β grows, stronger differences between both models arise. In the left column we have α/β = 1/5 and the dynamics of both systems are quite close (for the 3 first rows). In the right column α/β = 3/5 and stronger differences appear, specially for S(t). The most striking difference between both systems can be appreciated in the picture located at the last row, second column, which corresponds to a small perturbation of the case when initially half of the population is immunized and the other half is susceptible to the infection. In this case, the SIR system predicts no outbreak (it is a stable fixed point), while the modified SIR system does predict it. Obviously, this is due to the fact that in the modified SIR system we are assuming that the recovered population has been removed (death, quarantine, etc) and therefore does not interact (thus not contributing to the so-called 'herd immunity', which is of course not attainable in this model). All these considerations can also be deduced from the phase space representation in Figures 2 and 3 .

However, it is important to stress that the epidemiologically most relevant scenario, at least for a new epidemic like the COVID-19 one, in which there is no immunized individuals at the beginning (or they are very few ones), is given by the third row (cyan). So, we can conclude that in this scenario, specially when the ratio α/β is smaller, the SIR and modified SIR systems present similar features, with the modified SIR model always predicting a larger infection peak that the SIR.

Effective containment explains subexponential growth in recent confirmed COVID-19 cases in China. Science (80-. )

Simulating the spread of COVID-19 via spatially-resolved susceptible-exposed-infected-recovereddeceased (SEIRD) model with heterogeneous diffusion

Size and timescale of epidemics in the SIR framework

Asymptotic behavior in a deterministic epidemic model

The Mathematical Theory of Infectious Diseases

An integrable SIS model

A contribution to the mathematical theory of epidemics

Effects of information-dependent vaccination behavior on coronavirus outbreak: insights from a SIRI model

Estimation of COVID-19 dynamics 'on a back-of-envelope': Does the simplest SIR model provide quantitative parameters and predictions?

Application of Symmetry and Symmetry Analyses to Systems of First-Order Equations Arising from Mathematical Modelling in Epidemiology

Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates

Hamiltonian structure of compartmental epidemiological models

Accurate closed-form solution of the SIR epidemic model

Path Integral Approach to Uncertainties in SIR-type Systems

Models for the spread of universally fatal diseases

Poisson structure of dynamical systems with three degrees of freedom

A Modified SIR Model for the COVID-19 Contagion in Italy

This work has been partially supported by Ministerio de Ciencia e Innovación (Spain) under grants MTM2016-79639-P (AEI/FEDER, UE) and PID2019-106802GB-I00/AEI/10.13039/501100011033, and by Junta de Castilla y León (Spain) under grants BU229P18 and BU091G19.