Accurate closed-form solution of the SIR epidemic model.

An accurate closed-form solution is obtained to the SIR Epidemic Model through the use of Asymptotic Approximants (Barlow et al., 2017). The solution is created by analytically continuing the divergent power series solution such that it matches the long-time asymptotic behavior of the epidemic model. The utility of the analytical form is demonstrated through its application to the COVID-19 pandemic.

There are several problems of mathematical physics in which the only available analytic solution is a divergent and/or truncated power series expansion. Over the past decade, a new approach has evolved to overcome the convergence barrier in series solutions. An asymptotic approximant is a closed-form expression whose expansion in one region is exact up to a specified order and whose asymptotic equivalence in another region is enforced. The remarkable feature of asymptotic approximants is their ability to attain uniform accuracy not only in these two regions, but also at all points in-between, as demonstrated thus far for problems in thermodynamics, astrophysics, and fluid dynamics [1], [2], [3], [4], [5], [6], [7]. The current need to model and predict viral epidemics motivates us to extend the application of asymptotic approximants to the commonly used Susceptible–Infected–Recovered (SIR) model. This model is formulated as a system of nonlinear ordinary differential equations. Although no exact analytic solution has yet been found for the SIR model, a convergent series solution may be formulated via the homotopy analysis method [8]. Here, we provide an alternative and simple analytic approach. Interestingly, the SIR model shares the same asymptotic features as boundary layer flow over a moving flat plate, for which asymptotic approximants have already been applied [4]. The analytic nature of the asymptotic approximant derived in what follows is advantageous. Model parameters may be extracted for available COVID-19 data via a least squares (or equivalent) technique without the need for an embedded numerical scheme.

The SIR epidemic model considers the time-evolution of a susceptible population, , interacting with an infected population, , where is time. This model is expressed as [9] with constraints where , , , are non-negative constant parameters [9]. Once (1c) is solved, the recovered population is extracted as: Eq. (1a) can be thought of as a standard collision model in a 2-order chemical reaction, where species and “collide” to deplete the population of to create the species . In this interpretation, r is a rate constant, which in practice may be reduced by population behavior such as “social distancing”. In the case where =0 in (1b), the system (1c) indicates that for all time. For , then, the number of infected are reduced in time in accordance with (1b), and it is seen that the parameter determines the rate of recovery of infected individuals. The omission of a negative term in (1a) is an implicit model assumption that the recovered population is no longer susceptible to the disease.

We now manipulate the system (1c) into an equivalent first-order equation to simplify the analysis that follows. Eqs. (1a), (1b) are divided to obtain Subsequent integration of (3) with respect to and application of the constraints (1c) yields Eq. (4) is substituted into Eq. (1a) to obtain where From Eq. (1c), the constraint on is: System (5c) is equivalent to (1c) to solve for and, once solved, the solution for may be obtained using (4), which may be integrated to find from (2).

The series solution of (5c) is given by The result (6d) is obtained by applying Cauchy’s product rule [10] to expand and in (5c). The expansion of is obtained by first applying Cauchy’s product rule to the identity and evaluating like-terms to obtain a recursive expression for the coefficients of the expansion of , given by (6d). The expansion of is subsequently integrated term-by-term to obtain the expansion of , whose coefficients are given by (6c). Although the series solution given by (6d) is an analytic solution to (5c), it is only valid within its radius of convergence and is incapable of capturing the long-time behavior of . This motivates the construction of an approximant to analytically continue the series beyond this convergence barrier.

The long-time asymptotic behavior of the system (5c) is required to develop our asymptotic approximant, and so we proceed as follows. It has been proven in prior literature [11] that approaches a limiting value, , as , and this corresponds to according to (5c). The value of satisfies Eq. (4) with as [11]

We expand as as follows: Eq. (8) is substituted into (5c) and terms of O() are neglected to achieve the following linearized equation where In writing (9b), the definition of in (5b) has been employed. Additionally, to obtain (9b), Eq. (7) has been used which eliminates all O(1) terms in the linearized system. The solution to (9b) is where is an unknown constant that can only be determined via connection with short-time behavior. Consistent with the assumptions made, we find such that as . Thus the long-time asymptotic behavior of is given by Higher order corrections to the expansion (11) may be obtained by the method of dominant balance [12] as a series of more rapidly damped exponentials of the form where . This long-time asymptotic behavior of successive exponentials mimics that of the Sakiadis boundary layer problem describing flow along a moving plate in a stationary fluid [4]. It is natural, then, to apply the Sakiadis approximant [4] to capture this asymptotic behavior while retaining the behavior given by (6a). The Sakiadis approximant imposes the exponential form of the long-time asymptotic behavior (11) for all time; the coefficients of the exponentials are determined by matching their short-time expansion to the known power series developed about in the form of (6a). However, here we find that a reciprocal expression that achieves the same behavior (11) (through its binomial expansion) converges faster than the original Sakiadis approximant.

The assumed SIR approximant is given by where the ’s are obtained by taking the reciprocal of both sides of (12a), expanding each side about , and equating like-terms. The coefficients of the subsequent reciprocal expansion of the left-hand side (that of ) are given by (6d). After equating like-terms of this expansion with that of the reciprocal of the right-hand side of (12a), one arrives at the following linear system of equations to solve for the values as where (12b) is a Vandermonde matrix whose inversion is explicitly known [14]. The SIR approximant (12c) is thus a closed-form expression that, by construction, matches the correct behavior given by (11) and whose expansion about is exact to -order. A MATLAB code for computing the coefficients is available from the authors [15].

Fig. 1

Analytical and numerical solutions to the SIR model (5c) where the susceptible (), infected (), and recovered () populations are plotted versus time, all in arbritrary units. (a) As the number of terms is increased, the series solution, denoted (given by (6a), dashed lines), diverges and the approximant, denoted (given by (12c), solid lines), converges to the exact (numerical) solution (’s). (b) The converged asymptotic approximant for is used to obtain and (from Eqs. (2), (4), respectively). The model parameter values and initial conditions , , , and are taken from a test case used in Khan et al. [8] to validate the homotopy analysis method.

Fig. 1a provides a typical comparison of the -term series solution (6d) denoted by (dashed lines), the -term approximant (12c) denoted by (solid lines), and the numerical solution (indicated by symbols). Note that the series solution has a finite radius of convergence as evidenced by the poor agreement and divergence from the numerical solution at larger times, even as additional terms are included. By contrast, the approximant converges as additional terms are included. For , the approximant is visibly indistinguishable from the numerical solution (obtained by forward differencing) with a maximum relative error on the order of the numerical time-step (here 10−2) over the time range indicated. Increasing the number of terms beyond does improve accuracy up to a point, but also increases the likelihood of deficient approximants for which the denominator can be zero for certain time values and specific values of . In general, the lowest number of terms that yields the desired accuracy is chosen to avoid this behavior. The convergence of the approximant with increasing is a necessary condition for a valid approximant. For the problems of mathematical physics to which we have applied asymptotic approximants [1], [2], [3], [4], [5], [6], [7], we have observed that convergence of approximants implies excellent agreement with numerical results. There is as-of-yet no proof developed that guarantees this result, but this interesting behavior has been a property of all approximants developed thus far. In Fig. 1b, the converged () asymptotic approximant for is used to obtain and (from Eqs. (2), (4), respectively) and is compared with the numerical solution for these quantities.

In Fig. 2, the approximant is applied to a case examined in Khan et al. [8] to model the 1966 bubonic plague outbreak in Eyam, England. In Fig. 3, the approximant is applied to COVID-19 data for Japan [13]. An increased number of terms in the approximant is required to achieve the same relative errors in Fig. 1, Fig. 2, Fig. 3. For all cases examined, we observe that this trend correlates with the breadth of the initial plateau.

Fig. 2

Analytical and numerical solutions to the SIR model (5c) where , , and are in units of people and is in months. All other notation and labels are the same as in Fig. 1. The model parameter values and initial conditions , , , and are taken from estimates of the 1966 bubonic plague outbreak in Eyam, England examined in Khan et al. [8].

Fig. 3

Analytical and numerical solutions to the SIR model (5c) where , , and are in units of people and is in days. All other notation and labels are the same as in Fig. 1. The model parameter values   0.0164 and were obtained via a least-squares fit between the asymptotic approximant and Japan COVID-19 outbreak data [13] (’s), using initial conditions (from the first point in the data set [13]) and . Here is January 22, 2020 (see main text for interpretation of the COVID-19 data).

Note that the reported COVID-19 outbreak data [13] in Fig. 3 is originally provided in terms of confirmed cases and recovered individuals per day. The difference between these two quantities is used as an approximation to compare with the quantity of the SIR model. It is acknowledged that the actual COVID-19 data is influenced by transient effects not included in the SIR model such as the exposure lag-time; these effects are incorporated in more sophisticated models such as SEIR [11]. The approximation of from COVID-19 data is not restrictive in the current context, as our purpose is to show the efficacy of the closed form approximant rather than assess the validity of the SIR model.

In Fig. 3, a least squares fit of the asymptotic approximant to data is used to extract SIR parameters and based on data from the initial stages of the COVID-19 epidemic in Japan. To do so, (4) is used to relate analytically to the solution for (here, the approximant ); note that , used in the approximant, is affected by these parameters explicitly according to (7). The value of is not provided in the data set [13], and a least-squares algorithm is ineffective at determining an optimal value. Here, we choose the value of to be twice that of the maximum value of approximated from the data, as it captures a typical curve shape for seen in applications of the SIR model [11]. In regards to the sensitivity of fitting parameters to the choice for , a difference in leads to roughly a difference in and a difference in . The fit is made especially simple owing to the analytical form of the approximant that obviates the need to embed the numerical solution in such an algorithm. The population of recovered individuals, , is extracted from the solution for by direct integration in accordance with (2). Note that the predicted curve for in Fig. 3, obtained solely by fitting data for , is in good agreement with approximations from COVID-19 data for the recovered population, and serves as a check on the consistency of the data and algorithm.

It is evident from the results presented here that an asymptotic approximant can be used to provide accurate analytic solutions to the SIR model. Future work should focus on whether the asymptotic approximant technique can yield a closed form solution to more sophisticated epidemic models.

CRediT authorship contribution statement

Nathaniel S. Barlow: Conceptualization, Methodology, Software, Formal analysis, Visualization. Steven J. Weinstein: Conceptualization, Methodology, Formal analysis.