key: cord-0715140-dh7vs5cr
authors: Sadurní, E.; Luna-Acosta, G.
title: Exactly solvable SIR models, their extensions and their application to sensitive pandemic forecasting
date: 2021-02-01
journal: Nonlinear Dyn
DOI: 10.1007/s11071-021-06248-y
sha: 2c31b684c7398f79024050bedaeb1d8aee29792f
doc_id: 715140
cord_uid: dh7vs5cr

The classic SIR model of epidemic dynamics is solved completely by quadratures, including a time integral transform expanded in a series of incomplete gamma functions. The model is also generalized to arbitrary time-dependent infection rates and solved explicitly when the control parameter depends on the accumulated infections at time t. Numerical results are presented by way of comparison. Autonomous and non-autonomous generalizations of SIR for interacting regions are also considered, including non-separability for two or more interacting regions. A reduction of simple SIR models to one variable leads us to a generalized logistic model, Richards model, which we use to fit Mexico’s COVID-19 data up to day number 134. Forecasting scenarios resulting from various fittings are discussed. A critique to the applicability of these models to current pandemic outbreaks in terms of robustness is provided. Finally, we obtain the bifurcation diagram for a discretized version of Richards model, displaying period doubling bifurcation to chaos.

History Epidemiological models based on first-order, linear and multidimensional differential equations can be traced back to Daniel Bernoulli [1] . The famous susceptible-infected-removed (SIR) and SIRD (plus dead) systems were developed by Sir Ronald Ross nearly a hundred years ago; see [2] . The former entails the use of nonlinear terms in the model expressions [3] [4] [5] with the aim of describing interactions between S, I and R populations, together with the introduction of obvious critical points at vanishing populations of infected or susceptible individuals. Since the introduction of those models in the past century, the possibility of establishing control mechanisms based on predictions has been part of public health policies. The relevance of sufficiently simple contagion dynamics of this kind cannot be underestimated, particularly in recent times. There has been a continued interest in this subject as shown by [6] [7] [8] , but modern and ongoing research based on SIR models can be found in [9] [10] [11] . Motivation It may appear at first sight that the required equilibrium points for any sensible dynamical model of epidemics would impose nonlinearities which are insurmountable for an exact analysis. It is well known that some qualitative features of the SIR models can be extracted from the defining equations themselves, but the precise evolution of populations in terms of known functions necessitates explicit solutions. It is also desirable that the resulting expressions be amenable to numerical evaluations to avoid the use of robust numerical algorithms for multivariable first-order problems. For our purposes, it is important to know whether nonlinearities hamper full integrability. Integrability in the sense [12] that there exist a sufficient number of integration constants that allow separability of multivariable differential equations, which is more general than the existence of explicit solutions. The lack of this quality would relegate the study to numerically obtained flows and unpredictability when the system is in a chaotic regime. Fortunately, we will be able to show that the standard versions of SIR models allow separability and a reduction to quadratures. (In mechanical systems, these would be the integrals for periods or lapses.) Moreover, we shall see that the associated integral transforms allow evaluation under reasonable assumptions and approximations. The integrability of these models is another motivating factor to carry out future work in the direction of determining and designing appropriate policies to ameliorate the effects of an epidemic. Goals and results Separable systems fall into the category of integrable models. In this paper, we shall address the conditions for separability in the aforementioned SIR model. It is sometimes mentioned [2] that either, explicit solutions of SIR do not exist or they are unknown. We assert that this is not entirely correct: we shall see from the full integrability of the model that a last quadrature related to the lapse function or time interval can be found. The latter is given in terms of known functions, and it solves the problem explicitly in its entirety by reverse substitution into all previous expressions for population functions S, I and R. Incomplete gamma functions will appear in the results. The only caveat here is that such lapse or period integral cannot be put in terms of elementary functions and only an infinite series of them can be offered. However, their form is sufficiently amenable for numerical evaluation.

In the realm of non-autonomous systems and the application of policies in real time, we shall address a generalization of SIR with a sliding infection parametere.g. controllable by quarantine-and show that this problem is integrable as well if the total population is conserved and whenever the control parameter of infection depends explicitly on accumulated cases of infection at a given time. A generalization of SIR to interacting populations will be given in this work, and we will see that the number of separation constants does not match the number of degrees of freedom and then, only in this case, we shall have a lack of explicit solutions.

As an important application of our results, we shall address the problem of curve fitting in the context of COVID-19 acquired data in the case of Mexico, during the year 2020. A simple generalization of the logistic model, the Richards model, may capture the general behaviour of the accumulated infection curve, but errors in data acquisition may give rise to poor forecasting a few weeks ahead; the stability of the system in terms of data uncertainties will be studied numerically and analytically. Structure of the paper In Sect. 2 we present a selfcontained formulation of SIR and SIRD models, together with their variations in both solvable and nonseparable cases. In Sect. 2.1 explicit solutions are presented, including a lapse integral in terms of gamma functions and equilibrium conditions. Section 2.2 contains a treatment with variable infection rates. The case of two interacting regions is discussed in 2.3, and in Sect. 2.4 the general case of many interacting regions and total distributions is presented. Section 3 is devoted to an application of Richards model to the case of Mexico's COVID unknown-19 outbreak in the year 2020. In Sect. 3.1 the explicit solutions of the model are presented, together with an analysis of robustness under data fluctuations. Section 3.2 discusses the various forecasting scenarios resulting from curve fittings to COVID-19 in Mexico. In Sect. 3.3 we study the dynamics of a discretized version of Richards model, and in Sect. 4 we draw our conclusions.

The total population P(t) is a differentiable function that can be divided into susceptible individuals S(t), infected individuals I (t) and removed individuals R(t). It is assumed that I can infect others contained in S, but R is inactive. This removed population may contain cured and dead cases as R = C + D, where C and D are also considered inactive. It is also assumed that neither sources nor sinks are present in P, so this quantity is conserved; some generalizations relax this supposition. In this manner, C and D are not part of the dynamical system and they can always be recovered by a statistical treatment of death percentage D = q R, C = p R such that p + q = 1. This leaves us with SIR instead of SIRD. The velocity at which contagion occurs can be estimated in controlled (experimental) situations when a few members of the population come into contact. For two individuals sharing a common space of a given radius the initial values R = 0, S = 1, I = 1, P = 2 start a process for which at time t i we end up with S = 0 and I = 2. Because of proximity, the quantity t i is a function of population density. Since the new infected population is 2, S = −1 and we have S/ t i = −1/ t i ≡ −δ, so our new infection constant is reciprocal to the infection time. This situation can be scaled to more infected and more susceptible individuals; if a number S 0 shares the space with I 0 = 1, they will pass to the next category in the same time t i and if more infected individuals are active I 0 > 1 the probability of contagion increases proportionally. This establishes, in the limit of small time steps, the law dS dt = −δS I.

The constant δ so defined can be corrected in experiments involving larger populations. Similarly, the infected I suffer a change I which comes from two contributions: a depletion term due to removal (death or recovery) modulated by some coefficient α in a time t r which is proportional to I and a gain in a time t i controlled by a coefficient β which competes with the loss in S. We have dI dt = β S I − α I.

Finally, the removed population can be fed only from I , assuming there is no preventive cure to turn S into R (specifically C) directly. Therefore, the growth velocity of R is proportional to I with a constant γ also measured by the reciprocal time of the removal process:

Even in the case where the proportionality factors were unknown, the conservation of P = S+I + R forces δ = β, γ = α, as can be checked by adding the expressions (1), (2) and (3) 

where

is the infection-to-recovery time ratio. The system (4) is autonomous and it is exactly solvable, as we now show.

Let us consider κ constant in a fixed time interval. The interval could be taken as [0, ∞) for asymptotic analysis. The relations (4) imply the quotients

which constitute separable ordinary differential equations with solutions

To extract the time behaviour S(τ ), I (τ ), R(τ ) we eliminate I (R) in the third relation of (4) using the second relation in (7):

The last relation above is an ordinary differential equation for R alone and can be separated; thus, it can be solved by quadrature to yield the so-called lapse integral:

Integral transform for the lapse function From the growth relation d R/dτ = κ I , we know that if κ = 0 and I = 0, i.e. away from the trivial equilibrium point, it holds that d R/dτ = 0 for all times in [0, ∞). As a consequence, (9) defines a monotonic invertible function and R(τ ) can be recovered. Furthermore, the integral transform in (9) can be calculated in two regimes: κ < 1 and κ > 1. If κ < 1, then e (R 0 −X )/κ < 1 as here X > R 0 . A geometric expansion of the integrand is viable and obtains

where is the incomplete gamma function. At lowest order when κ → 0, one gets from either Eq. (9) or Eq. (10)

where we use κβ = α, c.f. Eq. (5) . This states that as t → ∞, R → P ≡ R max and all the population P is transferred to the removed state, as expected from a large infection rate β. On the other hand, if κ > 1, then ζ ≡ |e (R 0 −X )/κ −1| < 1 by expanding the exponential. In this case, a geometric series in ζ is plausible, and (9) leads once more to exponential integrals conforming a series of gamma functions:

To lowest order (e (R 0 −X )/κ ∼ 1), both Eqs. (8) and (9) yield

which implies R → R 0 + I 0 at infinity, a less drastic behaviour than the previous case (11) . To further distinguish between these cases, we note that κ also represents the susceptible population needed to reach the maximum infection, as can be inferred from the vanishing derivative in (4) . This means that κ 1 in realistic models, where infection peaks take place and S 1 simultaneously. Equilibria Before treating the case of jumping values of κ during a controlled outbreak, it is important to recall the critical behaviour of the involved quantities. From (4) we see that the infection peak occurs when S = κ, and if we use (7) the peak value I max results in

which is meaningful when I 0 +S 0 ≥ κ (log(S 0 /κ) + 1) (this is always the case if P is large). If κ = S 0 , then I max (min) = I 0 + S 0 − S 0 = I 0 and the peak occurs trivially at the initial condition. On the other hand, as κ increases, the function I max (κ) decreases when κ < S 0 ; indeed, if the susceptible number S decreases with the dynamics, S 0 is larger than the critical population S = κ and from (14) we find that the last two terms are negative, reducing the value of I max (κ) and explaining the famous curve flattening.

In order to recover the total deaths, we look now at R → R ∞ as τ → ∞. By direct time integration [see Eq. (4)] , we have

Equation (14) shows that I (τ ) depends on κ; hence, R ∞ is also κ-dependent. From the definition of κ = Eq. (5) and the meanings of α and β, modification of either β (by changing quarantine duration) or α (vaccine or cure) has an important effect on the total number of deaths D ∞ = q R ∞ . Remark 1 It can be shown, using (14) and (15) , that the artificial acceleration of infection (increase of β, thus decrease of κ) increases the final number of total deaths.

The quantity R ∞ can also be estimated from the final equilibrium point I ∞ = 0 by solving

A sequence of control parameters Suppose now that a new value κ is achieved at t = t f for all times in [t f , ∞). The previously computed quantities undergo a modification

A piecewise definition of κ consists of a set of constants {κ n } at intervals [t n , t n+1 ]. Then, new predictions are obtained for the dreadful I max :

Suppose now that I n = I (t n ) < I max (n) for a sufficiently large t n , i.e. at least one peak has already occurred. The condition for the emergence of a new peak of comparable proportions I max (n + 1) > I n is then

Using the total population P = S n + I n + R n = S n−1 + I n−1 + R n−1 in the relation above, the removed differential R = R n − R n−1 is shown to satisfy

The new parameter κ n is then given in terms of the old parameter and the death differential D/q. This quantity represents a comparison between consecutive regimes with infection peaks, and D > 0 implies worse scenarios in the evolution. In order to understand the effect of κ on D, let us define

From the dynamical equations, dS/dt < 0 regardless of κ, so we infer S n−1 > S n and we may bound D from below as

if f (x) < 0, and

x > e, so the infection is either too weak (large κ) or we are at the end of the outbreak (small S), which is not interesting for our analysis. In order to estimate D in the previous cases, it suffices now to study the sign of f in the inequalities above. We have x, x ∈ (0, e);

The first possibility yields 1 > κ n−1 /S n−1 > κ n /S n > κ n /S n−1 and so κ n < κ n−1 < S n−1 as expected, i.e. a smaller κ worsens the infection. However, the second possibility obtains eS n > κ n > S n κ n−1 /S n−1 , which represents a failed quarantine unless κ n grows beyond eS n , despite being bounded below by κ n−1 .

More realistic models of variable social behaviour can be studied. A non-autonomous system with κ = κ(τ ) can be proposed and solved explicitly in some cases. From the conservation of P, we are always left with a two-component system, as R is a redundant function of S and I :

As we saw before, it is always possible to re-parameterize using one of the SIR variables. For example, let our new time T be the monotonically increasing area under the curve I with time t as independent variable. For constant κ, this would be tantamount to parameterizing the dynamics with R alone. We have

, with

Then, (23) becomes

An iteration scheme can be applied to any nonautonomous model for general M(T ), similar to a Dyson series (physics [13] ) or a Volterra series (used by Levy, Wiener and others [14, 15] ). From (25) and the specific case at hand (26), we may find the explicit results S(T ), I (T ):

The remaining problem here is to find explicitly the function T (τ ) or its inverse. Nonlinear, second-order equation for lapse We expect T (τ ) to be governed by a non-autonomous equation, expressing that such kind of systems are not always separable, despite the fact that the existence of solutions is guaranteed by the Picard iterative method. The nonlinear, non-autonomous equation for T (τ ) can be obtained by differentiating with respect to τ the following expression

producing thus

This equation results from the application of the chain rule and the second relation in (27) . It is easy to see that for constant κ, (29) reduces to a separable equation whose quadrature coincides with (9) . Quadrature Let us solve (29) by separation of variables.

Here, it is important to note that κ(τ ) is assumed to be prescribed but, in contrast, κ(T ) needs the function T (τ ). We might think, however, that the application of some health policy rests on the knowledge of accumulated infections at time τ , which entails the use of κ(T ) instead of κ(τ ) as a specified control parameter. Then, by looking at (27) and (28) one has dτ = dT

and since κ(T ) is known by assumption, the quadrature

is a solution to (23) by direct substitution of the inverse T (τ ) into (27).

Example with an explicit form of the lapse Suppose a successful policy modelled by κ(T ) = κ 0 e T is applied. The lapse (31) is

where the values of κ 0 are such that the discriminant A 2 −4BC < 0. The function T (τ ) is recovered in terms of the functions tan and log. For a positive discriminant, tanh can be used instead. The maximum value of T is reached when τ → ∞, and it is determined by the condition A + Be −T + Ce T = 0. Comparisons In Fig. 1 we show the SIR evolution in weeks. The evolution of the dynamical variables was obtained by numerically integrating Eqs.(1)- (3) and their behaviour is consistent with our analytical expressions. A reasonable recovery time is 7 days, so 1/γ = t r defines our time scale (for COVID-19 it is two weeks [16] ). In compliance with [17] , numbers around 90-100 infected per 10 5 inhabitants are reported as world averages; these are used in our simulations: we take as initial conditions S 0 = 10 5 , I 0 = 10 2 , R 0 = 0 in all cases. We illustrate in panels (b) to (d) the fact that the asymptotic R ∞ , S ∞ are particularly sensitive to variations of κ. In fact, the area under the curve I is not constant with κ, as shown in (a) by plotting R.

In Fig. 2 we show the temporal evolution of S, R and I when κ also depends on time as a piecewise constant parameter and compare it with the temporal evolution of S, R and I when κ is constant throughout the whole epidemic. Recalling that κ = t i t r = α β , an increase in κ may be due to: (1) Increase in the removal rate α (i.e. a shorter time to recover) by, for example, improving the medical treatment at early stages or to a quicker time to die, due to perhaps some change in the environmental conditions, or (2) a decrease in the infection rate (β), by imposing a quarantine policy, for example. The first two weeks show a quickly evolving outbreak with κ −1 = 4 × 10 −3 , until a policy is applied, increasing κ by a factor of 4. From weeks 2-5 the infection curve seems under control, but a relaxed regime starting at week 5 is capable of generating a new peak in I . The revival occurs around week 7, and as a consequence, the asymptotics R ∞ increase D ∞ assuming constant q. Thus, the effect of increasing κ from week 2 to week 5 reduces notably the infection peak at the expense of a second milder regrowth and an increase in the duration of the epidemic. At the end the quarantine the removed population decreases by about a half, implying a drastic reduction in the number of deaths at the end of fifth week. However, right after week 5, the removed population picks up again and asymptotically it is only about 10 per cent smaller than the case of a constant κ throughout the whole outbreak. Figure 2 also shows that the effect of a quarantine depends on the initial and final date of its application.

When two regions come into contact by the exchange of susceptible and infected individuals, two species i = 1, 2 of quantities S i , I i , R i must be considered. It is reasonable to assume that infections between populations of separate regions are not possible, i.e. S 1 → I 2 and S 2 → I 1 are forbidden processes. However, S 1 ↔ S 2 and I 1 ↔ I 2 are plausible exchanges. There are two types of migratory trends that can be studied: the depletion of one population into the other and oscillatory dynamics between regions. It is important to note that in the limit of perfect isolation each set of SIR functions evolves as dictated by (4) . On the other hand, in the limit of negligible infectious rate, the system should decouple into independent migratory subsystems S 1 , S 2 , I 1 , I 2 and possibly R 1 , R 2 if removed patients are all recovered and allowed to travel. As before, the total population is conserved here, but the condition now reads P 1 + P 2 = P, with

The simplest way to take into account these considerations is by introducing an interaction term in our dynamical system

The functions F, G, J, K must comply with conservation of P and, in the absence of infection, the conservation of S = S 1 + S 2 , I = I 1 + I 2 and R = R 1 + R 2 . Therefore, F = −G and J = −K . The change of partial populations is thus d P 1 /dτ = F + J = −d P 2 /τ . Autonomous vs non-autonomous interaction Now we distinguish two important cases: Autonomous (timeindependent) interaction and non-autonomous (ad hoc) external signal. In the first case we have two subcases, namely a) logistic-type depletion and b) harmonic behaviour. For a) the interactions are

with λ and the intensity parameters given by the migration rate. The quantities σ and are critical points of the isolated migration model, associated with equilibrium populations. is the Heaviside step function, which prevents the undesired result S i , I i < 0 for any τ .

For the subcase b), a harmonic motion can be achieved by employing

where ω and are the frequencies and Sg is the signum function. That this is so, can be inferred from a simple change of variable in typical the harmonic systemẋ = ωy,ẏ = −ωx, where x 2 + y 2 is constant. Indeed, X = x 2 and Y = y 2 leads toẊ = 2ω √ |XY | = −Ẏ .

When the interactions are modelled by external "signals", a subcase b) can also be considered. Oscillations can be treated with

where a, A are oscillation amplitudes.

Remark 3 External signals such as (36) or damped signals with factors e − τ are useful for reproducing revivals. Bimodal curves for I i indicate cycles of infection parameterized by the frequency .

Interacting models such as (35) display rigidity for critical frequencies. A variety of numerical methods based on Runge-Kutta algorithms diverge at finite times.

The total population of two interacting regions does not behave as a full SIR, even in the absence of migration or exchange. Under the change of variables S = S 1 + S 2 ,

and the small subsystem is

Systems comprising two or more regions are, in general, non-separable.

Results For regions interacting with autonomous harmonic terms such as (35), we can show that the infection peak can display anomalous behaviour in different scenarios. The typical outbreak with initial conditions S 1 (0) = 10 5 = S 2 (0), I 1 (0) = 10 2 = I 2 (0), R 1 (0) = 0 = R 2 (0) is shown in Fig. 3 panel (a) . There we see a shoulder in the infection curve (black thick line) and an anomalous revival in the second susceptible population (blue thick line) whereas in panel (b) another set In such a case we also see an infection peak with a shoulder (black thin line) for the first population, but a drastically decreasing I 2 (t) (black, thick line). A more conspicuous oscillatory effect can be achieved by nonautonomous terms in (36): The results in Fig. 4 show an infection curve with notorious damped oscillations for parameters I 1 (0) = 2 × 10 3 = I 2 (0), S 1 (0) = 10 5 = S 2 (0), R 1 (0) = 0 = R 2 (0). The frequencies are ω = = 1/3 (in units of γ ). For the second type population, an infection curve with multiple "humps" in thick black is visible. Lastly, the results of a migratory depletion model in (34) are displayed in Fig. 5 . The difference between panels (a) and (b) comes from a larger migration rate for (b), albeit the use of similar initial conditions. As a result, the infection curve for population 1 (black thin line) has a larger peak in (b), in compliance with a migratory trend from region 2 to 1. 

Our previous remark establishes that separate regions produce a non-trivial joint effect even when they are isolated. Also, in the presence of migration, the evolution of total SIR populations of a country or large region unfolds differently as compared to non-interacting regions.

Uncorrelated case Let us suppose that during an outbreak a set of disconnected regions evolve under different, uncorrelated, conditions. The resulting individual curves will have similar shapes but different time shifts and peak intensities. For example, the infection processes could have different starting points, different initial conditions and varying values of κ. It is then natural to consider a full distribution made of individual populations with a weight factor W given by the local population fraction (e.g. France in the world: 

The integration has a smoothing effect in global measurable quantities. For instance, in a strict quarantine, the resulting curve for world populations in major scale events (see COVID-19 in [17] ) has a regular behaviour in limited time windows, as if modelled by a global SIR with effective κ. However, each region has a distinct behaviour governed by a local κ and entirely different sets of curves. With this, we would like to stress that the practice of performing a posteriori curve fitting does not necessarily validate the correctness of dynamical models locally.

Interacting case The set of equations for fully interacting regions become dS a dτ = −S a I a + 

Note that the second and third relations contain t and not τ in the limit. Once more we recall here that the absence of infection reduces to a model of circulation with conserved quantities SIR. More details can be introduced if a is replaced by coordinates r and F a,a , I a,a are replaced by generalized distributions containing ∇ r δ (3) − r ) , and so on, which give rise, upon integration, to local operators such as ∇, ∇ 2 , etc. From here we may obtain as special cases the diffusion equations and continuity conditions in the form of Fick's law.

The logistic model and its generalizations are useful in the description of growth for many applications, including biology. Phenomenological uses were originally proposed by Richards [18] , based on historical discussions by Gompertz [19] and previous work by Von Bertalanffy [20, 21] . Such logistic functions can be viewed as a particular reduction of SIR models to a single variable, i.e. the accumulated population A, which is the time integral of I under the condition P = S + I . This corresponds to the situation "once infected, always infected". The reduction gives rise to an integrable model by quadrature. Moreover, it allows to integrate the flow in terms of elementary functions for a variety of nonlinear systems that generalize the autonomous logistic equation. One of our aims here is to see how capable is the model, with constant parameters throughout the whole outbreak, to give reasonably good fits to the data and to express the model's sensitivity to parameter variations and data fluctuation; particular attention shall be paid to the total number D and its fluctuations. This can be easily done with analytic solutions at our disposal. Richards' model for the accumulated number of infected population A is defined by the flow

where the accumulated population A is given in terms of new infections I new

The parameters are r : modulation of the infection rate, K = A ∞ : final accumulated count of infected population , and a plays the role of recovery coefficient; the larger a is the faster the infection ends. The parameter a also determines the asymmetry of the curve of the daily new infections-e.g. a = 1 corresponds to a symmetric curve around A = K /2-with fixed points A = 0, A = K . The allowed values of the parameters are determined by the existence of such fixed points; we have r > 0, a > 0 and r < 0, −1 < a < 0, satisfying the condition ra > 0. In particular, a ∈ (−1, 0) generalizes the standard logistic map, where a = 1.

Assuming constancy of the parameters r, a, K , the autonomous system has solutions of the form

which translate into

Robustness, data fluctuations It is straightforward to compute the fluctuations in the predictions of final accumulated numbers of infected patients. The problem of interest is to analyse the error in forecasting due to bad data at a given time. Assuming that such statistical error at time τ is δ A, then the constant K changes to K + δ K . To see this, one simply solves for K in (45) as a function of A and τ and proceeds with the partial variation of A as independent variable, according to

Since A → K as τ → ∞, we expect δ K → δ A ∞ , as can be verified explicitly. Small errors in data produce then the variation

where the first fraction has been identified as the function K (A, A 0 ) obtained from (45). Expression (47) supports the concept that final predictions have different sensitivities at different times. For example, the variation on K due to a variation (bad data) of A at δτ very large (that is, near the end of the pandemic) is easily seen to be δ K = K δ A ∞ A ∞ = δ A ∞ , a negligible amount. In contrast, if the variation occurs at early times ( τ ≈ 0), then the impact on the prediction for K is δ K = δ A 0 K A 0 , which can be very large even if δ A 0 is small. Since this expression for δ K is a monotonously decreasing function of time, the errors due to bad data at some time t gradually become less important as t increases.

Richards model

Here, we try Richards model to fit Mexico's COVID-19 data [22] for forecasting and as an illustration of the sensitivity to small variations in data, obtained at various times. See also [23] . First, we obtained analytical expressions for relevant quantities, such as time of peak infections, number of daily new infections at its peak, etc. Integration of (42) with initial condition A 0 gives

where

Solution A(t) = (45) follows from (48). It is interesting to determine the time t p at which the speed of infection reaches its peak. This occurs when F(A; r, a, K ) reaches its maximum value, i.e. when ∂ F ∂ A = 0:

It follows that

Hence,

Note that the term (1 + a) −1/a in Eq. (50) is a monotonically increasing function of a, growing 0 at a = −1, reaching 1/e at a = 0 and approaching 1 at a → ∞ . The quantity (1 + a) −1/a can be roughly approximated by the straight line 1 e + e−2 2e a. So, A p grows roughly linear with a. Note that (50) limits the range of a to [−1, ∞) and, complementary, (52) indicates that negative values of a are allowed provided r is negative simultaneously with a. Also note that there is no problem with a being negative in the expression for t p = (52) since B 0 automatically becomes negative as a becomes negative It is also useful to determine the temporal evolution of daily new cases. In particular, we wish to know, for a given set of values of (a, r, K ) when and what is the highest number of daily new cases. This information can be obtained numerically using the solution (42) and computing A(t) = A(t) − A(t − 1). It can also be obtained by noting that F(A(t); a, r, K ) , where t = 1. This approximation is excellent since the function A(t) is very smooth for all t. For the particular time t = t p , the value F (A(t p ) ; a, r, K ) = r A p (1 − (A p /K ) a ) gives the maximum increase in daily new cases and can readily be obtained using (50):

The fitting parameters (a, r, K ) of solution (45) depend on the choice of precision of the fitting and on the number of data used to perform the fitting. Fitting precision is controlled by the maximum number of iterations M given as an input in a subroutine that incorporates the method of least squares. Increasing M improves the fitting up to a certain value M, after which the values of the parameters remain almost unchanged, according to a prescribed tolerance. Numerically we found M = 1150. Since the data are not very reliable due to various factors, such as the lack of actual laboratory testing in all suspected cases and the delay of gathering daily data, it is not indispensable to make as precise a fitting as possible. In this context, we may regard various different fitting precisions as possible different epidemic scenarios.

In what follows, we shall fit the solution (45) of Richards model to the data of cumulative number of infected people of COVID-19 in Mexico (data gathered from [17] www.worldometers.info/coronavirus). The first case of infection was detected on February 28, 2020, the day number 1. Table 1 shows the values of parameters (a, r, K ) and the predictions t p , A p and Ap, obtained from 8 representative fittings for two sets of data: one from day 1 till day 118 (June 24, 2020) and the second from day 1 till day 134 (July 10, 2020).

The first four fittings (I + , I − , II + , II − ) used the data of the first set and the last four cases (III + , III − , IV + , IV − ) used the second set of data. The subindex in each case refers to the sign of the parameter a (both signs are allowed, see discussion below (52)). Further, fitting I + differs from fitting II + in that the latter does not include the first 25 days of the data. The same distinction is true for the pairs (II + , IV + ). This was done in order to Here and elsewhere n i (n f ) is the initial (final) day of the fitting data test if predictions could be improved by resting importance to the initial data and to test how much the fitting parameters change with time. Comparing the values of K , A p , A, t p and t f for the pair (I + , II + ) shows that, indeed, disregarding the initial 25 data points yields somewhat larger values of the number of total accumulated infections and the duration of the epidemic. As we shall see below, the larger values of K , t p , A p makes the extrapolation of the fitting II + come closer to the data values for a few days after the end of the fitting, the day 118. The root mean square deviation σ , a measure of how much the fit deviates from the data, is appreciable smaller for the case II + than for the case I + , indicating that case II + is a better fitting. This improvement on the fitting correlates with an improvement on the short term prediction. The same occurs between cases III + and IV + ; this time with 134 as the final day of fitting; here, σ changes from 1856 to 1562 as the first 25 days are ignored.

In Fig. 6 we plot A(t) = (50) for the 8 cases listed in Table 1 together with the actual data. Figure 7 is a zoom of Fig. 6 . The above-mentioned cases correspond to the four lowest curves, I + (blue solid line), II + (black solid line), III + (red solid line) and IV + (magenta solid line). These illustrate that even though cutting off the first 25 day of fitting does improve the predictions, it does it only for a week or so after the fitting final date.

It is evident that all 8 fittings give similar results for the first stages of the pandemic but their asymptotic behaviour is quite different, fanning out into two sets. Specifically, the set that diverts very soon from the data corresponds to those fittings with positive values of a and r . The second set of fittings that have all negative values of a (simultaneously, negative values of r ) is able to predict the data correctly and for a much longer time.

The reason negative values yield better predictions is that the data for the daily new cases A are far from being symmetric about its peak value A p . It is symmetric when the peak A p is K /2. Richards model allows for the asymmetry through the parameter a. As Eq. (50) shows, the flow F(A; r, α, K ) = (42) is symmetric only for a = 1 ( the standard logistic flow). As a becomes less than 1, the peak moves to the left of A = K /2 and as a → 0 + , A p → K /e. So, the solution (45) would be unable to provide good fittings with a > 0 if the data demands that A p be smaller or equal to K /e. This argument is consistent with the ratios computed for all the fittings we have done. In particular, the ratios A p /K for the cases I + , II + , III + and IV + are, respectively, 0.36793, 0.36793, 0.36791, and 3.6794; these are barely larger than 1/e = 0.36787 . . .. The other four cases have negative a (and r ), and the ratios are little smaller than 1/e, namely 0.3512, 0.34317, 0.35791 and 0.34313 for I − , II − , III − and IV − , respectively. Note that as a becomes negative, r becomes negative also in order to preserve the concavity of the flow function (42)

In conclusion, Figs. 6 and 7 show that the fittings with negative values of a and r render not only a better fitting but their extrapolations after the final fitting date are able to predict quite well the data for about 2 months. Still, there is no certainty that the predictions will continue to agree with data beyond a few days after the day 188 ( September 2, 2020) but one may 

For any given country or region the generation of data for epidemics is updated daily, so it is interesting to see what differences may occur for the discretized version of Richards model. Rescaling the cumulative number of infections A(t) to the total accumulated number of infections K : A → x ≡ A/K the discretized version of (42) is the mapping

where we have taken the constant lapse time t between gathering of data as our unit of time. We can further simplify the flow (42), by rescaling the time by the infection coefficient r , t → τ ≡ r × t, leading to

Note that we have changed the notation for the scaled number of infections x = I /K to y = I /K since x(t) is in general not equal to x(r × t) ≡ y(τ ). Its corresponding mapping is

Here, our unit of lapse time is τ = r t. We now list the basic properties of this simpler homeomorphism Based on these properties, we conclude that the interesting interval of α is (0, α t ) since in this interval all orbits are bounded and the period one fixed point y p = 1 loses its stability at α = 2. What remains to see is what kind of dynamics occurs from α = 2 to α t . Figures 8 and 9 succinctly display the rich dynamics of the mapping H α (y) in the interval 1.8 ≤ α ≤ α b . The interval α ∈ (0, 1.8) is not included since all orbits converge to the fixed point of period one y p = 1, as mentioned in item vi above. This bifurcation diagram shows not only the (predicted) instability of the period-one fixed point for α > 2 but that the dynamics undergo the period doubling route to chaos. However, it should not be surprising given that our simple mapping belongs to the universal class of unimodal mappings with negative Schwarzian derivative in the interval [0, y e ]. [24, 25] .

Remarks As regards the implementation of this simple map to the epidemics, one may also fit the data to the orbit (cascade) y n = H α (y n ). Of course, we would want a scenario where the orbits converge to one single outcome, which implies that α should be less than 2. Small values of α, say α = 0.2 yield sigmoid-like evolutions, Fig. 10 shows the evolution of y(n) and daily new cases y n = y(n + 1) − y(n) for two representative small values of α = 0.2 (blue lines and dots and α = 0.2 (black lines and dot) with same initial conditions, x(1) = 10 −6 .

The two parameter model F(x, α, r )= Eq. (42). A similar analysis can be carried out for the two parameter mapping (42). It can be shown that the maximum occurs in three regions in parameter space, namely region 1) r > 0, α < −1, region 2) r < 0, −1 < α < 0 and region 3) r > 0, α > 0. Calculations of the Schwarzian derivative of F(x, α, r ) show that it is negative in large regions where the critical value of F(x, α, r ) is a maximum, so it is expected that it will also display the period-doubling route to chaos. The detailed analysis of this model shall be reported elsewhere.

In this paper, we have provided integrability conditions and exact solutions for a sufficiently simple model of SIR. In small steps, we improved the model so as to encompass the variation of parameters with time and the interaction between regions, including migration Note that the plots for y n were multiplied by a factor of 10 for better visibility in a single graph of susceptible and infected individuals. A system with more than two regions is non-separable in any of its non-trivial versions; the exceptions, of course, are the absence of population exchange or the absence of infection. A variety of effects were discussed for infection revival under migratory conditions and variable quarantine policies. With respect to applicability, some comments are in order. In particular, Richards model (a generalized logistic model) was employed to fit Mexico's COVID-19 data up to a certain date of the epidemic. As a result, a variety of scenarios were obtained, consistent with strong fluctuations in the final number of casualties. It was possible to find a family of curves fitting the data up to a certain date reasonably well, but with strong variations in the final outcome, months ahead. In this respect, we have seen how small variations in parameters produce large deviations in projected number of accumulated cases, and therefore, using D = q R, also in the number of casualties. This exhibits the challenge of reliable metrics for SIR and logistic type systems; a great deal of work and an impressive display of modelling methods in real time has been seen during the year 2020. To name a few important examples, we refer the reader to [26] [27] [28] . For more recent studies, see [29] [30] [31] . However, while acknowledging the merit of these valuable works, we should underscore the fact that forecasting and nowcasting are strongly challenged by the careful determination of phenomenological parameters, especially the reproduction number R 0 . Small inaccuracies may lead to very large errors, not because of chaoticity (which is absent) but because of exponential growth and the subsequent exponential depletion.

In this respect, the existence of exact solutions allows for a systematic study of robustness, understood as error propagation under model variation and uncertainties. Interestingly, in the branch of quantum dynamics and complex modelling applied to physical systems, the concept of fidelity is sometimes mentioned as a useful measure, as it consists of an autocorrelation function between solutions after evolution takes place under two slightly different dynamical systems. We believe that a similar analysis can be carried out in population dynamics with the results provided in the present work.

Finally, it is important to stress the role of inverse problems, as opposed to direct problems, typically addressed in this setting. The former consists in the determination of a model from a desired behaviour, whereas the latter deals with the determination of the evolution from a given model. Inverse problems tap into exactly solvable systems and explicit formulas because of invertibility, a key concept guaranteed by the implicit function theorem.

Author contributions Both authors conceived the project, wrote the manuscript and revised mathematical and numerical calculations.

Funding There was no funding specific for this work.

Daniel Bernoulli's epidemiological model revisited

The SIR model and the foundations of public health

Three basic epidemiological models

The mathematics of infectious diseases

The basic epidemiology models: models expressions for r0, parameter estimation and applications

Discussion: the Kermack-McKendrick epidemic threshold theorem

Analysis of sir epidemic models with nonlinear incidence rate and treatment

Dynamics of an epidemic model with quadratic treatment

COVID-19 pandemic in India: a mathematical model study

Application of mathematical modeling in public health decision making pertaining to control of COVID-19 pandemic in India

COVID-19 Healthcare demand and mortality in Sweden in response to non-pharmaceutical mitigation and suppression scenarios

Modern Quantum Mechanics, Revised edn

Bibliography on Volterra Series, Hermite Functional Expansions, and Related Subjects

Nonlinear System Theory: The Volterra-Wiener Approach

WHO Director-General's opening remarks at the media briefing on COVID

Reported cases and deaths by country, territory, or conveyance

A flexible growth function for empirical use

On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies

Stoffwechseltypen und Wachstumstypen

Quantitative laws in metabolism and growth

Report on the forecast of cumulative COVID-19 cases in Mexico

Iterated Maps of the Interval as Dynamical Systems

Philip Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcatios of Vector Fields

Early dynamics of transmission and control of covid-19: a mathematical modelling study

Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in Wuhan, China: a modelling study

The effect of control strategies to reduce social mixing on outcomes of the covid-19 epidemic in Wuhan, China: a modelling study

Analytical features of the SIR model and their applications to covid-19

Inversion of a SIRbased model: a critical analysis about the application to covid-19 epidemic

A SIR model assumption for the spread of covid-19 in different communities

The authors declare that they have no conflict of interest.