Exact and approximate analytic solutions in the SIR epidemic model
EXACT AND APPROXIMATE ANALYTIC SOLUTIONS IN THE SIR EPIDEMIC MODEL
Mário Berberan-Santos Institute for Bioengineering and Biosciences Instituto Superior Técnico Universidade de Lisboa 1049-001 Lisboa Portugal [email protected]
Abstract
In this work, some new exact and approximate analytical solutions are obtained for the SIR epidemic model, which is formulated in terms of dimensionless variables and parameters, reducing the number of independent parameters from 4 ( I , S , , ) to 2 ( i = I / S and R = S / ). The susceptibles population is in this way explicitly related to the infectives population using the Lambert W function (both the principal and the secondary branches). A simple and accurate relation for the fraction of the population that does not catch the disease is also obtained. The explicit time dependences of the susceptibles, infectives and removed populations, as well as that of the epidemic curve are also modelled with good accuracy for any value of R using simple functions that are modified solutions of the R → limiting case (logistic curve). It is also shown that for small i ( i < 10 -2 ) the effect of a change in this parameter on the population evolution curves amounts to a time shift, their shape and relative position being unaffected. Keywords:
SIR epidemic model, Kermack-McKendrick model, epidemic dynamics, Lambert W function.
Highlights • The susceptibles population is explicitly related to the infectives population using the Lambert W function (both principal and secondary branches). • The time dependence of the susceptibles, infectives and removed populations, as well as that of the epidemic curve are described with good accuracy using simple functions. • A change in the initial number of infectives on the population evolution curves is shown to amount to a time shift (provided the initial number of infectives is much smaller than that of the susceptibles).
1. Introduction
The propagation of new viral infectious diseases can be described, to a first approximation, according to the SIR model developed by Kermack and McKendrick [1-6]. In this deterministic compartmental model, three classes of individuals are considered: Susceptibles, S , who can catch the disease; Infectives (infected and infectious), I , who can transmit the disease to susceptibles; and Removed, R , who are former infectives that either recovered from the disease and became immune to it (permanently or in the time frame of the analysis) or died because of the disease. It is usually assumed that no members of the Removed class are initially present. Given the timescale of a typical epidemic, birth, death (by other causes) and other demographic processes can usually be neglected. It is further assumed that the latency period is negligible, implying that an individual, once infected, immediately becomes infective. The model also assumes a perfect mixing of Susceptibles and Infectives, as a result of free motion (extensive interconnectedness) within the closed system. With these assumptions, the model can be written as the following pair of evolution steps: S I I + ⎯⎯→ (1) I R ⎯⎯→ (2) The first step accounts for the infection process that occurs via close contact of an infective with a susceptible, and is the infection rate constant. It reflects the intrinsic frequency of contacts and the probability of transmission upon contact: according to the theory of diffusion-influenced collisional processes, it can be written as = pc , where c is the encounter rate for full diffusion control and p is the probability of transmission upon encounter [6]. The second step describes the immunization or death of infectives that thus become members of the Removed class, and is the removal rate constant, representing the probability of transition I → R per unit time. The nonlinear system of differential equations corresponding to the SIR mechanism is [1]: dS S Idt = − , (3) dI SS I I Idt = − = − , (4) dR Idt = , (5) where S , I and R are population densities (number of class elements per unit area). As R does not appear in the first two equations, the problem is effectively reduced to a system of two (autonomous) differential equations. The dynamics is determined by four parameters: (i) Initial number of susceptibles S ; (ii) Initial number of infectives I (assumed to be suddenly mixed with the susceptibles at time zero). I is usually much smaller than S (it can be as low as a single individual in the entire system); (iii) Rate constant that reflects the likeliness of disease transmission upon contact between an infective and a susceptible; (iv) Rate constant that reflects the average time needed for an infective to be removed. The inverse of is the average infectious period [2-5], which is also the average duration of an infective. Another quantity of interest is the epidemic curve, C [2-5] C S I = , (6) which is the number of new infectives per unit time, as follows from equation (4). The purpose of this work is to obtain new explicit relations for the SIR model (formulated in terms of dimensionless variables) both exact and approximate, including simple analytic solutions for the time-dependence that are valid for any R . The paper is organized as follows: In Section 2, the SIR model is formulated in terms of dimensionless variables and the susceptibles population is explicitly related to the infectives population using both the Lambert W function and an approximate yet accurate function (sum of exponentials). The maximum value of the epidemic curve is also explicitly obtained. In Section 3, the explicit time dependence of the susceptibles, eq. (38), and infectives, eq. (41), as well as that of the epidemic curve, eq. (45), are obtained in terms of approximate analytic functions whose respective parameters are obtained in Section 4 as a function of R . The introduced functions also allow to define analytically the respective maximum times. In Section 5, the effect of i ( i << 1) on the time-dependence is shown to amount to a time shift. Finally, in Section 6 the main results are summarized.
2. Results
Upon application of Laplace transforms, it follows from equations (3) and (4) that exp exp t tI I C = − + − , (7) ( ) R I I I C S I S I = = − + = + − + , (8) where stands for the convolution between two functions. The epidemic curve thus acts as the impulse function that defines the infectives time evolution. The system of differential equations can be rewritten in terms of five reduced (dimensionless) quantities: Ss S = , (9) Ii S = , (10) Rr S = , (11)
00 0
SR S = = , (12) t = , (13) where R is the basic reproduction number of the infection [2-6]. This quintet is preferable to another set, not entirely dimensionless, previously proposed [7]. Equations (3)-(5) become: ds R s id = − , (14) ( ) di R s id = − , (15) dr id = . (16) The initial values of s and i are 1 and i , respectively. Usually, i << 1. The reduced epidemic curve, c , is c R s i = , (17) In this way, the evolution of the system is fully defined by a single initial condition, i , and by a single dynamic parameter, R . This parameter has the meaning of the average number of susceptibles infected by a single infective at the beginning of the epidemic [2-7]. Values up to 15 are known [8], while for COVID-19 in pre-pandemic conditions it was between 2.4 and 2.5 [9,10]. The R value reflects not only the intrinsic aspects of the disease, but also the population specific conditions (frequency of contacts, probability of transmission). The threshold (or critical) number of susceptibles, is c s R = . (18) If R < 1 ( s c > 1), the number of infectives monotonically decreases with time (cf. equation (15)) and the epidemic cannot proceed (‘herd immunity’ situation) [1-6]. On the other hand, if R > 1 ( s c < 1), the epidemic ensues, with the number of infectives increasing until s c is reached. Owing to the nonlinear terms, no explicit analytical solution exists for the system of differential equations (14)-(16), and the time evolution must be obtained numerically or using approximate equations. A few analytical results are nevertheless possible. Elimination of time from equations (14) and (15) yields dids R s = − , (19) whose integration gives the (phase space) relation between i and s [1-6]:
11 ln i i s sR = + − + , (20) depicted in Fig. 1.
Fig. 1
The phase space plot: i vs. s according to equation (20). A value of i = 0.1 was used. The number next to each curve is the respective R . The dashed line ( R = ) corresponds to the logistic curve situation. The arrow of time is also shown. The maximum value of i , i c , occurs for s = s c = R , as follows from equation (20), and takes the value (see Figure 1) ( )
11 1 ln c i i RR = + − + . (21) It is useful to solve equation (20) explicitly for s . The result is (see Appendix A): ( ) (1 )1 00 i i R c s W R e s sR − + −− = − − , (22) ( ) (1 )0 00 i i R c s W R e s sR − + − = − − , (23) where W (x ) is the Lambert function computed for the principal branch [11] and W -1 (x ) is the Lambert function computed for the secondary branch [11]. Using equations (22) and (23), the value of S can be found directly for different points of the epidemic, instead of solving numerically equation (20). In particular, the final value of s (value at the end of epidemic), s , follows from equation (23) by setting i = 0, ( ) (1 )0 00 i R s W R eR − + = − − , (24) as previously obtained from the final size equation [12,13]. Assuming that i << 1, as is usually the case, equation (24) becomes ( ) R s W R eR − = − − , (25) For R higher but very close to 1, equation (25) reduces to 2/ R -1. For large R , it becomes exp(- R ). Instead of equation (24), a simple and accurate formula (average relative error < 1% if i < 10 -3 ) and with the correct asymptoptic behaviour can be used, ( ) ( ) ( ) exp 2.462exp 1.851 8.798exp 3.580 ( 1.1) s R R R R = − + − + − . (26) The dependence of s on R is plotted in Figure 2 (see also Figure 1). As discussed above, there is no epidemic for R < 1 and the number of susceptibles does not decrease significantly ( i << 1). The fraction of the population that does not catch the disease is about 40% for R = 1.5, 20% for R = 2, and 6% for R = 3. For R > 5 virtually the whole (initially susceptible) population is infected. Fig. 2
The fraction of susceptibles not catching the disease, s vs. the basic reproduction number R , as given by equations (25) and (26). It is assumed that i << 1. Also shown is the critical value s c , equation (18) (dashed line). Another important parameter that can be obtained is the maximum value of the epidemic curve, c * (Appendix B):
1* * * c R s s R = − , (27) where s * stands for the number of susceptibles when the production of new infectives is at the maximum, ‘epidemic peak’ (Appendix B), ( ) [1 (1 ) ]1 00
1* 22 i R s W R eR − + +− = − − , (28) where W -1 ( x ) is the Lambert W function computed for the secondary branch. For large R equation (27) becomes ( i << 1): * 4 Rc = . (29) 0 Finally, the areas under both i and c curves can be obtained from equations (16) and (14), respectively, and are found to be nearly identical ( i << 1): ( ) 1 i d i s = + − . (30) ( ) 1 c d s = − . (31) In this way, for large R ( R > 5, say) and very small i , both i and c curves are area normalised. The respective peak values, equations (21) and (29), respectively, agree with the fact that for large R the curve i ( ) attains a stable shape (exponential decay function), with a peak value of 1+ i , whereas the curve c ( ) is increasingly narrow and approaches a delta function.
3. Simple functions describing the time-dependence
Obtention of simple closed-form solutions for the time-dependence of the quantities s , i , r and c has thus far eluded all efforts, unless R is very close to 1 [1-3]. Recently, Barlow and Weinstein [14] presented an accurate closed-form solution using asymptotic approximants. However, the number of terms needed can be quite large and no explicit relations can be obtained for characteristic quantities such as the maximum times. A different approach is presented here, starting from the limiting results that correspond to R → . In fact, for large R , one has R s >> 1 until s is very close to 0. In this way, virtually all S is transformed into I before R starts to form. Equation (15) can therefore be approximated by di R sid = , (32) implying that i + s = i + s in this time range. The solution of the system of equations (14) and (31) is: 1 ( ) is i Ri += + + . (33) The respective i ( ) is given by the logistic equation [15], ( )
00 0 0 0
1( ) exp 1 ii i i R i += − + + . (34) Finally, c ( ) is ( )( ) exp 11( ) 1 1 exp 1 i Rc R i i Ri + = + + + , (35) and peaks at ( )
00 0 ln* 1 ii R = − + , (36) with the value ( )
20 0 * 14
Rc i = + , (37) compare equation (29). Equation (36) shows the effect of the initial fraction of infectives on the dynamics. The lower their number, the higher the induction time. Full time range approximate solutions for the SIR model and for any R are now obtained using eq. (33) as the starting point for a trial function representing s ( ). The trial function is written as ( ) ( )
11 exp as s s a b += + − + , (38) 2 where parameters a and b should reduce to 1/ i and to (1+ i ) R , respectively, in the limit of large R . The parameter s is given by equation (24). The only modification in the mathematical form is thus to allow for s > 0. A simple functional form is kept, allowing to analytically compute integrals and characteristic parameters such as those for the maxima of s and c . All these quantities can indeed be obtained in closed form from the trial function equation (38). Equation (15) gives
00 0 ln ( ) 1 i R s u dui = − , (39) and therefore ( ) b i aR si a b a a e + = − + + + , (40) or ( ) b a ei i ea e − += + , (41) with
1( 1)(1 ) a s ba s + −= + − , (42) Rs a b = − + . (43) This function has the correct long-time behaviour for large R , decaying as exp(- ). When R → , parameter → b , whereas parameter →
1, and the decay of i becomes exponential for all times. The maximum value of the infectives curve is attained at 3 ( )
11 ln 11 c R ab R s −= + − . (44) Note that the alternative computation of i ( ) from equation (14) is not effective, as a and b are no longer independent and the smoothing effect on deviations by means of integration is nonexistent. The epidemic curve is obtained from equations (1), (17) and (38), ( ) ( ) ( )( ) exp1 1 exp bdsc s a bdt a b = − = − + + , (45) and the time at which the maximum (‘epidemic peak’) is attained is ln* ab = , (46) hence the respective peak value is ( ) ( ) ( ) a b sc b s aa + −= − , (47) compare equations (29), (36) and (37). In practice, c ( ) is best computed using the defining equation (17). Finally, the r curve can in principle be computed either using equation (8) or from ln sr R = − . (48) This last relation, obtained by the elimination of time from equations (14) and (16) [3], is less convenient as it amplifies the error of s when this quantity is close to zero.
4. Determination of parameters a and b Curve fitting for the determination of parameters a and b (fixed i and R ) can be performed in two ways. In the first method, the s and i curves are obtained by numerical integration, and the fitting is carried out in the time domain in a global way (i.e., minimizing a single sum of squares of deviations and using the same parameters for both curves). In the second method, fitting is carried out in the phase space, no numerical integration being necessary. In this case the exact i ( s ) curve is given by equation (20), and the fitting i ( s ) curve is given by equation (41), the reduced time being computed as a function of s with equation (38), rewritten as ( ) sab s s −= + + − . (49) The dependence of parameters a and b with R for i = 0.001 are shown in Figures 3 and 4 (the numerical values are given in the Supplementary Information). Parameter a tends to 1/ i when R → , as discussed in connection with equation (38). On the other hand, it is observed that it tends to 1 when R →
1. Parameter b follows an almost linear dependence with R , becoming close to 0 for R = 1, as expected. A similar pattern is observed for other values of i (up to 0.01; no higher values tested). Fig. 3
Plot of parameter a as a function of R , obtained from the phase space fitting of i vs s ( i = 0.001). The initial value is a (1) = 1, and the limiting value ( R → ) is 1/ i = 1000, see equation (38). Fig. 4
Plot of parameter b as a function of R obtained from the phase space fitting of i vs s ( i = 0.001). It is observed that the dependence is almost linear and that b (1) is close to zero. Fitting for several values of R (1< R < 1000) confirms they allow a very good representation of the results. Representative results of phase space fittings are shown in Figures 5-8. Fig. 5
Phase space plot for R = 1.3 and i = 0.001 using equation (20) and the respective fitting using equations (41) and (49). Fitted parameters: a = 74.40 and b = 0.2782 ( = 0.6481 and = 2.017). Fig. 6
Phase space plot for R = 3.0 and i = 0.001 using equation (20) and the respective fitting using equations (41) and (49). Fitted parameters: a = 274.5 and b = 1.683 ( = 1.789 and =1.683). Fig. 7
Phase space plot for R = 6.0 and i = 0.001 using equation (20) and the respective fitting using equations (41) and (49). Fitted parameters: a = 522.8 and b = 4.577 ( = 4.588 and =1.310). Fig. 8
Phase space plot for R = 10.0 and i = 0.001 using equation (20) and the respective fitting using equations (41) and (49). Fitted parameters: a = 702.7 and b = 8.593 ( = 8.593 and =1.165). Minor deviations for some values of s exist for intermediate values of R , as seen in Figures 6 and 7. Nevertheless, the overall shapes and the values and position of the maxima are quite well reproduced in all cases. It is possible to use more complicated functions to get better fits for s (and therefore for the remaining quantities), but at the expense of introducing additional parameters and of losing analytical power, for instance ( ) ( ) ( )( ) a ps s s a b += + − + , (50) with an additional parameter ( p ), gives very good fits for s ( ), almost eliminating the somewhat too fast decay observed after the middle point, see e.g. Figures 10 and 11. The chosen set of equations represents a compromise between accuracy and simplicity, allowing to obtain explicit expressions for the peak values and respective times, for instance equations (44), (46) and (47). 8 Fig. 9
Plot of the s , i and c curves as a function of time for R = 1.3 and i = 0.001, both obtained by numerical integration (color) and by phase space fitting (black, dashed). * = 15.4 and c = Fig. 10
Plot of the s , i and c curves as a function of time for R = 3.0 and i = 0.001, both obtained by numerical integration (color) and by phase space fitting (black, dashed). * = 3.34 and c = Fig. 11
Plot of the s , i , r and c curves as a function of time for R = 6.0 and i = 0.001, both obtained by numerical integration (color) and by phase space fitting (black, dashed). * = 1.37 and c = Fig. 12
Plot of the s , i , r and c curves as a function of time for R = 10.0 and i = 0.001, both obtained by numerical integration (color) and by phase space fitting (black, dashed). * = 0.763 and c = Figure 12 clearly demonstrates the impulse-response relationship between c ( ) and i ( ), as expressed by equation (7). In particular, after c ( ) (essentially) ceases, i ( ) decays exponentially.
5. Change in i : time shift effect The numerical results presented in the previous section refer to i = 0.001. From fittings with other values of i it is concluded that parameter a has an inverse dependence with i , for constant R , whereas parameter b practically does not change with i ( i < 0.01), especially for R > 2. However, for small i ( i < 0.01), there is no need to evaluate a and b parameters as a function of i . The effect of changing this parameter is equivalent to a time shift of the curves, without altering their shape or relative position. Indeed, equation (39) can be rewritten as 1 ( ) exp ( ) 1 i i R s u du = − , (51) hence, ( ) exp ( ) 1 exp ( ) 1exp ( ) 1 exp ( ) 1 . i i R s u du R s u dui R s u du R s u du − = − − = = − + − (52) Maclaurin series expansion of s gives ( ) 1 ... s R i = − + , (53) and therefore, for << 1/( R i ), one may use s ( ) = s (0) = 1, and equation (52) becomes ( ) ( ) exp 1 exp ( ) 1exp ( ) 1 i i R R s u dui R s u du −− = − + − = = + − , (54) where the new initial value, '0 i , is ( ) '0 0 0 exp 1 i i R = − . (55) Comparison of eq. (54) with eq. (51), rewritten as ''0 00 '( ') exp '( ) 1 i i R s u du = − , (56) 2 gives '( ) ( ) s s = + , (57) and therefore, the time shift relations are obeyed both by s and i , '( ') ( ) s s = , (58) '( ') ( ) i i = , (59) where ' = − , and the time shift is '00 0 iR i = − . (60) If the a and b parameters for i = 0.001 are used, then a time shift is needed for other values of i . For instance, if i = 10 -7 and R = 3, then the shift (time required to reach ' 30 i − = ) is = 4.6, implying an increase of the induction time.
6. Conclusions
The SIR epidemic model was formulated in terms of dimensionless variables and parameters, thus reducing the number parameters from four ( S , I , , ) to two ( i , R ). The susceptibles population was explicitly related to the infectives population using the Lambert W function (both the principal and the secondary branches), eqs. (22) and (23). A simple and accurate relation for the fraction of the population that does not catch the disease was obtained, eq. (26). The maximum of the epidemic curve (epidemic peak) was also obtained, eq. (27). The explicit time dependences of the susceptibles, infectives (and removed) populations were modelled with good accuracy for any value of R using simple functions matching the limiting case R → (logistic equation), eqs. (38) and (41). From these, the epidemic curve and the peak time are derived, eqs. (45) and (46). It was also shown that for small i ( i < 10 -2 ) the effect of a change in this parameter on the population 3 evolution curves amounts to a time shift, eq. (56), their shape and relative position being unaffected. Acknowledgements
This work was suported by Fundação para a Ciência e a Tecnologia (FCT), Portugal [project UIDB/04565/2020].
Competing interests statement
The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 4
Appendix A: Derivation of equations (22) and (23)
Equation (20) can be rearranged to give ( ) ln 1 s R s i i R − = − + − , (A1) or ( ) exp 1 R s s e i i R − = − + − . (A2) Defining y R s = − ( ) exp 1 y y e R i i R = − − + − , (A3) hence, from the definition of the Lambert W function [11], ( ) ( ) exp 1 y W R i i R = − − + − , (A4) and finally ( ) ( ) s W R i i RR = − − − + − . (A5) Given that the argument of the function is negative, the appropriate branch of W (principal or secondary) must be defined. Let ( ) exp 1 x R i i R = − − + − . (A6) Then ( ) s W xR = − . (A7) 5 This function is plotted in Figure A1 for R = 3 and i = 0.1 . It is seen that the secondary branch, W -1 ( x ), applies for t [0, t max ] and the principal branch W ( x ) for t [ t max , + [. Fig. A1
Variable s vs x , for R = 3 and i = 0.1, showing the two branches of W ( x ). s max represents the value of s when i attains its maximum value, i max . Appendix B: Derivation of equations (27) and (28)
The rate of production of new infectives ( epidemic curve ) is - ds/d , and the maximum rate occurs when d sd = . (B1) It follows from equations (B1), (14) and (15) that i * = s * - = s * - s c , (B2) where s * and i * stand for the number of susceptibles and infectives when the production of new infectives is at the maximum ( c *). The maximum rate is thus
1* * * c R s s R = − . (B3) Using equations (B2) and (20) it is obtained that
00 0 s s iR R − − + + = . (B4) It follows from equation (B2) that s * > s c . The solution of equation (B4) is then ( ) [1 (1 ) ]1 00
1* 22 i R s W R eR − + +− = − − , (B5) where W -1 ( x ) is the Lambert W function computed for the secondary branch (see Appendix A). The dependence of s* with R ( i << 1) is plotted in Figure B1. The reduced parameter starts from s * = 1 for R = 1, but quickly approaches the asymptotic value 1/2 as R increases and is reasonably constant for R > 5. 7 Fig. B1
The dependence of s* with R ( i << 1) . A simple formula, with an accuracy better than 0.1%, and with the correct asymptoptic behaviour is ( i << 1): ( ) ( )
1* 3.40 exp 2.17 0.143exp 0.247 2 s R R = − + − + . (B6) Given that s* approaches 1/2 for large R , it follows from equation (B3) that c* also tends to a constant value for high R ( i << 1), * 4 Rc = . (B7) Indeed, for large R the logistic curve applies in this time range and s and i cross at s * = i * = 1/2. 8 Literature cited [1] Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proc R Soc Lond A 1927; 115: 700-721. https://doi.org/10.1098/rspa.1927.0118 [2] Murray JD. Mathematical Biology. 2 nd Supplementary information Recovered a and b parameters and related quantities as a function of R ( i = 0.001) R a b s c -5 -5 -6 -6 -7 -7 -7 -8 -8 -9 -9
50 918.0009 48.47956 1.93×10 -22 -44 -66 -87 -218498.748 1.003524 0.026287 0.013829 1000 994.139 999.177 0 999.177 1.00183 0.013821 0.006908