Analytic solution of the SEIR epidemic model via asymptotic approximant
Steven J. Weinstein, Morgan S. Holland, Kelly E. Rogers, Nathaniel S. Barlow
aa r X i v : . [ q - b i o . P E ] J un Analytic solution of the SEIR epidemic model via asymptotic approximant
Steven J. Weinstein,
1, 2
Morgan S. Holland, Kelly E. Rogers, and Nathaniel S. Barlow a) School of Mathematical Sciences, Rochester Institute of Technology, Rochester, NY 14623,USA Department of Chemical Engineering, Rochester Institute of Technology, Rochester, NY 14623,USA (Dated: 1 July 2020)
An analytic solution is obtained to the SEIR Epidemic Model. The solution is created by constructing asingle second-order nonlinear differential equation in ln S and analytically continuing its divergent powerseries solution such that it matches the correct long-time exponential damping of the epidemic model. Thisis achieved through an asymptotic approximant (Barlow et. al, 2017, Q. Jl Mech. Appl. Math, 70 (1), 21-48)in the form of a modified symmetric Pad´e approximant that incorporates this damping. The utility of theanalytical form is demonstrated through its application to the COVID-19 pandemic.Asymptotic approximants have been successful at pro-viding analytical solutions to many problems in mathe-matical physics . Like the well-known Pad´e approxi-mant , they are constructed to match a primary seriesexpansion in a given region up to any specified order. Un-like Pad´e approximants, however, the form of an asymp-totic approximant is not limited to a ratio of polynomi-als, and its structure is chosen to enforce the asymptoticequivalence in a region away from the primary series ex-pansion. By increasing the number of terms in an asymp-totic approximant, it converges to the exact solution inthese two regions – as well as at all points in between.Convergence is certainly a necessary condition for a validapproximant; although there is yet no proof, convergentapproximants match the numerical solutions of systemsexamined thus far .The COVID-19 outbreak motivates the application ofasymptotic approximants to epidemiology models. Themethod has seen recent success in providing a closed-formsolution to the Susceptible–Infected–Recovered (SIR)model . Here, we extend the method to the commonlyused Susceptible–Exposed–Infected–Recovered (SEIR)model. This model is formulated as a system of non-linear ordinary differential equations, for which no exactanalytic solution has yet been found. The analytic natureof the asymptotic approximant derived in what follows isadvantageous, in that the accuracy and computationalexpense are not affected by the duration of the epidemicprediction; the form is built such that it is accurate in t ∈ [0 , ∞ ) and can be evaluated at any specific time with-out the need for numerical marching. Depending on theduration, it may be beneficial to replace a numerical so-lution with the approximant within a fitting algorithmthat extracts SEIR parameters. En route to the approx-imant, we also present an alternative formulation of theSEIR model as a single 2nd-order nonlinear differentialequation in ln S . This form enables an efficient series so-lution about t = 0, asymptotic expansion as t → ∞ , andmay itself prove attractive for future analysis. a) Electronic mail: [email protected]
The SEIR epidemic model considers the time-evolutionof a susceptible population, S ( t ), interacting with an ex-posed population, E ( t ), and infected population, I ( t ),where t is time. This model is expressed as dSdt = − βSI (1a) dEdt = βSI − αE (1b) dIdt = αE − γI, (1c)with a removed population (recovered + deaths), R ( t ),evolved by dRdt = γI (1d)and constraints S = S , E = E , I = I , R = R at t = 0 . (1e)In (1), β , α , γ , S , E , I , and R are non-negativeconstant parameters . Along with initial conditionsfrom (1e), the solution for S , E , and I may be first ob-tained from (1a) through (1c) and the solution for R subsequently extracted using (1d).We now manipulate the system (1) into an equivalent2nd-order equation in ln S to simplify the analysis thatfollows. Equations (1a) and (1b) are added to obtain dSdt + dEdt = − αE. (2)Solving (1c) for E and substituting into (2) then leads to d Idt + ( γ + α ) dIdt + α dSdt + αγI = 0 . (3)(1a) is rewritten as I = − β d ln Sdt (4)and substituted into (3) to arrive at the 3rd-order equa-tion d ln Sdt + ( γ + α ) d ln Sdt − αβ dSdt + αγ d ln Sdt = 0 . (5)Equation (5) may be integrated to yield d ln Sdt + ( γ + α ) d ln Sdt − αβS + αγ ln S = C, (6)where the integration constant C = αγ ln( S ) − αβ ( E + I + S ) (7a)is obtained by evaluating the left-hand side of (6) at t = 0using (1c), (1e). and (4). The form of (6) suggests thatthe variable substitution f = ln S be made, and the resultis d fdt + ( γ + α ) dfdt − αβe f + αγf = C (7b)where, from (1e) and (4), f = ln S , dfdt = − βI at t = 0 . (7c)Once (7) is solved for f , S is extracted as: S = e f . (8a)The solution for I follows directly from (4) and (8a) as I = − β dfdt . (8b)After substituting (8b) into (1d), integrating, and apply-ing the constraint (7c), R is expressed as: R = R − γβ ( f − ln S ) . (8c)Lastly, the conservation of S + E + I + R provides asolution for E as E = E + I + S + R − I − S − R, (8d)as seen by adding equations (1a) through (1d), integrat-ing in t , and applying (1e).The series solution of (7) is given by f = ∞ X n =0 a n t n , a = ln S , a = − βI (9a) a = [ C − ( α + γ ) a + αβS − αγa ] / a n +2 = αβ ˜ a n − ( γ + α )( n + 1) a n +1 − αγa n ( n + 2)( n + 1) , n > a n> = 1 n + 1 n X j =0 ( n − j + 1) a n − j +1 ˜ a j , ˜ a = S . (9d)The result (9) is obtained by the standard procedure ofinserting (9a) into (7) and finding a recursion for the co-efficients by equating like-terms. It is thus necessary toobtain the expansion of the nonlinear term e f ≡ S in (7).To do so, we solve for the coefficients of S = P ˜ a n t n byapplying Cauchy’s product rule to the chain-rule result f ′ S = S ′ and evaluating like-terms; this leads to therecursive expression given by (9d). Although the seriessolution given by (9) is an analytic solution to (7), it isonly valid within its radius of convergence and is inca-pable of capturing the long-time behavior of the system.This motivates the use of an approximant to analyticallycontinue the series beyond this radius.The long-time asymptotic behavior of the system (7)is required to develop our asymptotic approximant, andso we proceed as follows. It has been proven in priorliterature that S approaches a limiting value, S ∞ , as t → ∞ , and this corresponds to I → f approaches a limiting value, f ∞ ≡ ln S ∞ , as t → ∞ . The value of f ∞ satisfies the following equation e f ∞ − γβ ( f ∞ − ln S ) − E − I − S = 0 (10a)in the interval f ∞ ∈ ( −∞ , ln γ/β ) . (10b)We expand f as t → ∞ as follows: f ∼ f ∞ + g ( t ) where g → t → ∞ . (11)(11) is substituted into (7b) (with (7a)), e g is replacedwith its power series expansion, and terms of O ( g ) areneglected to achieve the following linearized equation d gdt + ( γ + α ) dgdt + (cid:0) αγ − αβe f ∞ (cid:1) g = 0 . (12)The general solution to (12) is g = ǫ e λ t + ǫ e λ t (13a) λ , = 12 (cid:20) − α − γ ± q ( γ − α ) + 4 αβe f ∞ (cid:21) (13b)where ǫ and ǫ are unknown constants and λ < λ < e f ∞ < γ/β from (10b). Thus the long-time asymp-totic behavior of f is given by f ∼ f ∞ + ǫ e λ t , t → ∞ . (14)Higher order corrections to the expansion (14) may beobtained by the method of dominant balance as a se-ries of more rapidly damped exponentials. However, thepattern by which the corrections are asymptotically or-dered is not as straightforward as that of the SIR model,provided in Barlow and Weinstein . In that work, anasymptotic approximant is constructed as a series of ex-ponentials that exactly mimics the long-time expansion.In the SEIR model, complications in the higher-orderasymptotic behavior arise from the competition betweenthe two exponentials in (13a). Here, we enforce theleading-order t → ∞ behavior given by (14) and makea more traditional choice for matching with the t = 0 ex-pansion (9). We create an approximant with an embed-ded rational function with equal-order numerator and de-nominator (i.e., a symmetric Pad´e approximant ), suchthat it approaches the unknown constant ǫ in (14) as t → ∞ , while converging to the intermediate behavior atshorter times. The assumed SEIR approximant is givenby f A,N = f ∞ + e λ t N/ X n =0 A n t n N/ X n =1 B n t n , N even (15)where the A n and B n coefficients are obtained such thatthe Taylor expansion of (15) about t = 0 is exactly (9).Note that, although a rational function is being usedin (15), it is not a Pad´e approximant itself. Pad´e approx-imants are only capable of capturing t n behavior in thelong-time limit, where n is an integer. The pre-factor e λ t is required to make (15) an asymptotic approximant forthe SEIR model. However, we may still make use of fastPad´e coefficient solvers by recasting (15) as a Pad´eapproximant for the series that results from the Cauchyproduct between the expansions of e − λ t and f − f ∞ ,expressed as N X n =0 n X j =0 ( − λ ) j j ! ˜ a n − j t n = N/ X n =0 A n t n N/ X n =1 B n t n , (16)where ˜ a = a − f ∞ and ˜ a n> = a n> . A MATLAB codeto compute the A n and B n coefficients of (15) (for given α , γ , β , S , E , I ) is available from the authors .The SEIR approximant (15) is thus an analytic expres-sion that, by construction, matches the correct t → ∞ behavior given by (14) and whose expansion about t = 0is exact to N th-order. A comparison between the approx-imant solution (15) and the numerical solution to (1) isprovided in figures 1-4 with the relative error for all fourcases provided in figure 5. The indicated error in figure 5is calculated by comparing S ( t ) to its accurate numericalsolution (assumed to be exact); curves showing the sameorder of accuracy are obtained when the other dependentvariables of the model are examined.Figure 1a provides a typical comparison of the N -termseries solution (9) denoted by f S,N (dashed lines), the N -term approximant (15) denoted by f A,N (solid lines), and the numerical solution ( • ’s). Note that the seriessolution has a finite radius of convergence as evidencedby the poor agreement and divergence from the numer-ical solution at larger times, even as additional termsare included. By contrast, the approximant convergesas additional terms are included. For N = 18, the ap-proximant is visibly indistinguishable from the numericalsolution on the scale of figure 1a. Figure 5a provides therelative error of the approximant for the data shown infigure 1. Increasing the number of terms beyond N = 18does improve accuracy up to a point, but a minimum er-ror barrier is eventually reached of O (10 − ) at N = 26;note that, to make this assessment, we take the maximumrelative error with respect to time for each N (the max-ima in figure 5a). For larger values of N , the maximumerror increases, and the approximant begins to diverge,i.e. there is an optimal value of N at which to truncatethe approximant. Asymptotic approximants can exhibitan optimal truncation as is often observed with asymp-totic expansions in general . We emphasize here that anumerical solution is not needed to assess convergenceof approximants to within their optimal truncation; con-vergence in the Cauchy sense (i.e., the distance betweenapproximants decreases with increasing N ) may be ex-amined. In addition to this issue, deficient approximantsare possible with increasing N due to zeros that can arisein the denominator of (15). Such approximants are ig-nored in assessing convergence. To avoid this behavior,the lowest number of terms that yields the desired accu-racy should be chosen. The convergence of the approxi-mant with increasing N (up until its optimal truncation)is a necessary condition for a valid approximant. In fig-ure 1b, the converged ( N = 18) asymptotic approximantfor f is used to obtain analytic solutions for S , E , I , and R from (8), which are compared with the numerical so-lution for these quantities. The approximant for N = 18agrees with numerics within the visible scale of the plot,with errors quantified by figure 5a.Figure 1 results described above correspond to a caseexamined in Rachah and Torres to model an Ebolaoutbreak. In figures 2, 3, and 4, the approximant is ap-plied to COVID-19 data for Yunan (China), Sweden,and Japan, respectively. Figure 5b-d provide the rel-ative error for these cases; the largest indicated valueof N in each figure (corresponding to dashed curves) isthe optimal truncation as discussed above for figure 5a.Note that we extensively surveyed the available COVID-19 data , and the results in figures 2-5 are representa-tive of the fits and variability in the number of termsneeded for convergence of the approximant up to its op-timal truncation.Note that the reported COVID-19 outbreak data isprovided in terms of confirmed cases, recovered individ-uals, and deaths per day. We use recovered + deaths asan approximation to the removed population R and use confirmed − recovered − deaths as an approximation to I in the SEIR model. It is acknowledged that the ac-tual COVID-19 data is influenced by effects not included FIG. 1. Analytical and numerical solutions to the SEIR model (1), where the susceptible ( S ), exposed ( E ), infected ( I ), andrecovered ( R ) populations are represented as a fraction of the total population and t is in units of days. (a) Solution shown interms of f ≡ ln S . As the number of terms N is increased, the series solution, denoted by f S,N (given by (9), dashed curves),diverges and the approximant, denoted by f A,N (given by (15), solid curves), converges to the exact (numerical) solution ( • ’s).Corresponding relative errors are provided in figure 5a. (b) The converged asymptotic approximant for f is used to obtain S , E , I , and R from (8) shown by solid curves and compared with the numerical solution (closed symbols). The model parametersvalues and initial conditions α = 0 . β = 0 . γ = 0 . S = 0 . E = 0 . I = 0 .
05, and R = 0 are taken fromestimates of Ebola virus propagation examined in Rachah and Torres . FIG. 2. Analytical and numerical solutions to the SEIR model (1), where S , E , I , R are in units of people and t is in days.All other notation and labels are the same as in figure 1, except R now also includes deaths. Corresponding relative errors areprovided in figure 5b. SEIR model parameters values and unknown initial conditions are obtained via a least-squares fit tothe Yunan, China COVID-19 outbreak data (open symbols). Best fit parameters are α =0.395031, β =0.00333, γ =0.0553093, S =142, and E =0. The initial conditions I = 44 and R =0 are taken directly from the data set at a chosen t = 0 (hereJanuary 28, 2020). FIG. 3. Analytical and numerical solutions to the SEIR model (1), where S , E , I , R are in units of people and t is in days.All other notation and labels are the same as in figure 1, except R now also includes deaths. Corresponding relative errorsare provided in figure 5c. SEIR model parameters values and unknown initial conditions are obtained via a least-squares fit tothe Sweden COVID-19 outbreak data (open symbols). Best fit parameters are α =0.041281, β =1.513332 × − , γ =0.004407, S =50306, and E =10015. The initial conditions I = 1743 and R =20 are taken directly from the data set at a chosen t = 0(here March 21, 2020). FIG. 4. Analytical and numerical solutions to the SEIR model (1), where S , E , I , R are in units of people and t is in days.All other notation and labels are the same as in figure 1, except R now also includes deaths. Corresponding relative errorsare provided in figure 5c. SEIR model parameters values and unknown initial conditions are obtained via a least-squares fit tothe Japan COVID-19 outbreak data (open symbols). Best fit parameters are α =0.2332207, β =2.040015 × − , γ =0.034334, S =15442, and E =0. The initial conditions I = 1649 and R =529 are taken directly from the data set at a chosen t = 0(here April 1, 2020). -14 -12 -10 -8 -6 -4 -2 -14 -12 -10 -8 -6 -4 -2 -14 -12 -10 -8 -6 -4 -2 -14 -12 -10 -8 -6 -4 -2 FIG. 5. Relative error of the approximant (15) for increasing N as a function of t (in days). The exact solution is taken tobe the numerical solution of (1), computed using the 4th-order Runge–Kutta scheme with a time-step of 10 − . The subfigures(a)-(d) correspond to the cases presented in figure 1-4, respectively. For all figures, N is taken up until optimal truncation isachieved, indicated by a dashed curve. The cusps in the figures have no physical meaning and simply indicate where the signof ( S A,N − S exact ) changes. in the SEIR model, and this can affect the ability of themodel to closely fit actual COVID-19 data. The data ap-proximations made here are to enable comparisons withmodel predictions. The ability of the approximant tomatch numerical results is unaffected by such approxi-mations. Disagreement between the model and epidemicdata after fitting is attributed to the applicability of theSEIR model and not the approximant.In figures 2-4, a least squares fit to I and R data isused to extract SEIR parameters α , β , γ and initial con-ditions S and E . To do so, the initial values of I and R are taken directly from the COVID-19 data set . Ad- ditionally, the time t = 0 is chosen such that disease hasprogressed to a point where initial trends are observed,so that curve shapes are consistent with those reason-ably predicted by the SEIR model. Adjustments suchas this have been well described in fits done in previ-ous work . The initial guesses for the iterative least-squares fit are taken from data fits for earlier times thanexamined here .Our results demonstrate that an asymptotic approxi-mant can be used to provide accurate analytic solutionsto the SEIR model. Future work should examine theability of the asymptotic approximant technique to yieldclosed-form solutions for even more sophisticated epi-demic models, as well as their endemic counterparts . N. S. Barlow, A. J. Schultz, S. J. Weinstein, and D. A. Kofke,J. Chem. Phys. , 204102 (2012). N. S. Barlow, A. J. Schultz, S. J. Weinstein, and D. A. Kofke,AIChE J. , 3336 (2014). N. S. Barlow, A. J. Schultz, S. J. Weinstein, and D. A. Kofke,J. Chem. Phys. , 071103:1 (2015). N. S. Barlow, C. R. Stanton, N. Hill, S. J. Weinstein, and A. G.Cio, Q. J. Mech. Appl. Math. , 21 (2017). N. S. Barlow, S. J. Weinstein, and J. A. Faber, Class. Quant.Grav. , 1 (2017). R. J. Beachley, M. Mistysyn, J. A. Faber, S. J. Weinstein, andN. S. Barlow, Class. Quant. Grav. , 1 (2018). E. R. Belden, Z. A. Dickman, S. J. Weinstein, A. D. Archibee,E. Burroughs, and N. S. Barlow, Q. J. Mech. Appl. Math. ,36 (2020). N. S. Barlow and S. J. Weinstein, Physica D , 1 (2020). G. A. Baker Jr. and J. L. Gammel, J. Math. Anal. Appl. , 21(1961). C. M. Bender and S. A. Orszag,
Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Pertur-bation Theory (McGraw-Hill, 1978). W. O. Kermack and A. G. McKendrick, Proc. Roy. Soc. LondonA , 700 (1927). H. W. Hethcote, in
Mathematical Understanding of InfectiousDisease Dynamics , edited by S. Ma and Y. Xia (World ScientificPublishing, 2008). A. Rachah and D. F. M. Torres, Math Method Appl. Sci. ,6155 (2017). “https://github.com/chebfun/chebfun/blob/master/padeapprox.m,”. P. Gonnet, S. G¨uttel, and L. N. Trefethen, SIAM Rev. , 101(2013). John Hopkins University CSSE, “Novel coronavirus (covid-19)cases,” https://github.com/CSSEGISandData/COVID-19. M. Peirlinck, K. Linka, F. S. Costabal, and E. Kuhl, Biomech.Model. Mechan. doi.org/10.1007/s10237-020-01332-5 , 1(2020). K. Linka, M. Peirlinck, F. S. Costabal, and E. Kuhl, Comput.Method Biomec. doi.org/10.1080/10255842.2020.1759560doi.org/10.1080/10255842.2020.1759560