Analytical solution of the SIR-model for the temporal evolution of epidemics. Part A: time-independent reproduction factor

Several recent studies [1–4] have demonstrated that the normal or Gaussian distribution function for the temporal evolution of the daily number of new cases (deaths, or alternatively infections) at time t due to the COVID-19 pandemic disease provides quantitatively correct descriptions for the monitored rates in many different countries during a single wave. If applied early enough at the beginning of a pandemic wave, within the regime of non-exponential growth, the Gauss model (GM) makes realistic and reliable predictions for the future evolution of the first wave. It has been argued that the assumption of a Gaussian time evolution is well justified by the central limit theorem of statistics [3], an agent-based model [4], a Taylor expansion [4], or as a special case of the SIR (susceptible-infectious-recovered/removed) model [2]. A motivation of the present manuscript was to provide more rigorous arguments in favor of using the GM to estimate the characteristics (peak time, height, and width) of a first epidemic wave well ahead of its climax.

To describe the time evolution of pandemic deceases many models have been developed which are based on the original deterministic compartmental SIR-model pioneered by Mc-Kendrick and Kermack [5] and its later variants such as the SEIR-model (for recent review see Estrada [6]). Important modifications to these mean-field rate equations include stochastic aspects to account properly both, for strong number fluctuations, driving the continuous phase transition of the system near the epidemic threshold, and also for spatially correlated clusters and spreading fronts caused by the disease transmission through nearest-contact infection [7, 8]. It is hoped that with these important modifications the so improved SIR models ultimately provide quantitatively correct descriptions for the forecast of pandemics in many countries, even beyond the exponential growth phase of the pandemic. The inclusion of the stochastic effects mentioned above is often done by combining numerical computer simulations with the mean-field rate equations of the SIR- or SEIR-models with constant values of the infection and reaction rates. Therefore these studies will also profit from more general analytical solutions to the SIR-model equations for time-dependent infection and recovery rates which we develop here.

On our way we discovered a new analytical solution of the standard SIR-model [5, 9–13] without vital dynamics describing the temporal evolution of the COVID-19 pandemic disease that applies to the whole range of parameters. This solution allows for an arbitrary time dependence of the infection and recovery rates but assumes that the ratio of the two rates is independent of time. We emphasize that the assumption of a constant value for the ratio k of the recovery to the infection rate is necessary to deduce an analytical solution for the SIR model equations. There is no reason why the ratio of these two rates should be constant, as the recovery rate μ(t) is largely determined by the nature of the disease which is not expected to change strongly over time, whereas on the other hand, the infection rate a(t) is affected by many factors, not the least mitigation measures introduced by the government. Even if we ignore the latter non-pharmaceutical interventions (NPIs) it is not obvious that the infection rate varies proportional to the recovery rate. Nevertheless we are convinced that our calculations based on a constant value of the ratio k are meaningful and important, basically for the following four reasons: first, our case generalizes earlier treatments where time-independent infection and recovery rates were adopted. Secondly, as our solution allows for an arbitrary time dependence of the infection and recovery rates, despite their constant ratio, it is possible for the first time to model analytically the influence of mitigation measures on the time evolution of epidemics as done in reference [14]. Thirdly, our analytical solution can serve as a benchmark for the verification of numerical solutions of the SIR-model and their variants. Fourthly, if the intrinsic time scales of the pandemic event and the taken mitigation measures are rapid enough we may use our analytical results derived for a constant value of k to study slowly time-varying ('adiabatic') time-dependences k(t) by inserting these slow time variations into our final analytical expressions for k. Such a Wentzel–Kramers–Brillouin-approximation has proven useful in many branches of propagation and transport theory of physics including quantum mechanics.

Moreover, in contrast to earlier work we will also calculate analytically the daily differential rate of newly infected persons resulting from the SIR-model which is the key quantity to compare with the monitored data in different countries, as the size of the currently infected, and not yet recovered, compartment is usually not known. Its asymptotic behavior, peak time and peak amplitudes will all be obtained analytically and exactly.

Besides providing analytic expressions for the quantities characterizing the solution of the SIR model, we derive a simple, accurate approximant that can be used in practice, and shares all relevant features with the exact solution, as we will show. This is a significant improvement compared with approaches, where the solution of the SIR model was for example expanded into a divergent but asymptotic series [15, 16], or where it had been obtained assuming inequalities that do hold only within a very limited range of SIR parameters, as we will show.

As a side-observation we find that the SIR model exhibits an asymptotic behavior at large times that is qualitatively different from the logistic model, while the difference between the two models still decreases with increasing reproduction factor. Because the SIR model is sometimes used with arbitrary initial conditions, we recall here that initial conditions must be interrelated if the SIR model is assumed to hold at all times. This is the scenario to be investigated here, for reasons to be discussed in detail. The present investigation can be extended to the many variants of the SIR model [17–25].

As the dynamics of the COVID-19 pandemic is much faster than the dynamics of births and deaths, the neglect of these demographic factors is well justified. The SIR system is the simplest of the compartmental models used for the mathematical modeling of infectious diseases. The considered population of N ≫ 1 persons is assigned to the three compartments S (susceptible), I (infectious), or R (recovered/removed). Persons from the population may progress between these compartments.

2.1. Basic equations

In a fixed population of individuals let S(t), I(t) and R(t) denote the susceptible, infected and recovered/removed fractions of persons involved in the infection at time t, so that

and because S, I, R are fractions, they must all reside within the interval [0,1]. If a(t) and μ(t) denote the semipositive time-dependent infection and recovery rates, respectively, the SIR-model is defined with the two dynamical equations [5, 9, 26]

where the dot here and in the following denotes a derivative with respect to t. The equation for the dynamics of R(t) follows from the sum constraint (1), i.e.

where we inserted equations (2) and (3). Equation (3) can be written as

so that

implying that equation (4) becomes

with the potentially time-dependent inverse reproduction factor

The solution and its analytic approximant to be derived in this work holds for arbitrary a(t), as long as μ(t) and a(t) remain proportional to each other. In that case k is a constant, usually denoted as inverse basic reproduction number, k = 1/r₀.

2.2. Initial conditions: all-time versus semi-time SIR

There are two qualitative very different cases to be considered that can be regarded as equally valid approaches. We refer to these cases as the all-time case (I) and the semi-time case (II).

For the all-time case (I), treated by Kendall [9] and part A of this work, the ratio k = μ(t)/a(t) is regarded as identically constant at all times, from t = −∞ to t = ∞. This implies, that one has to use boundary conditions S(−∞) = 1, I(−∞) = 0, R(−∞) = 0 as the epidemic must not have existed at t = −∞. For the special case of k = 0, μ(t) = 0 and R(t) = 0 at all times according to equation (7), while equations (1)–(3) reduce to I(t) = 1 − S(t) and the simple logistic differential equation determining all fractions analytically in terms of a(t) (appendix F).

For constant k > 0, equation (7) implies the important relation [9]

valid for any choice of t and t'. Upon inserting the stated boundary conditions at t' = −∞ into equation (9), S(t) and R(t) are simply related via S(t) = e^−R(t)/k or R(t) = −k  ln  S(t). A special case of this relationship is

In view of equation (1), there remains the freedom to just specify one single initial condition at a certain t = 0, where data may have become available. One may regard the time t = 0 when the existence of the pandemic wave in the society is realized and monitoring of newly infected persons starts. Besides k there is thus a second parameter of the SIR model, which we denote by the positive ɛ as it usually represents a small number, but our results are valid for any ɛ. We define ɛ via the susceptible fraction of the population at time t = 0,

There is no freedom for the remaining R(0) and I(0) if k is known. Inserting the above boundary conditions into equation (10), and making use of equation (1) one has

involving a function in f_k (ɛ) that is going to frequently occur in this work

Here R(0) and I(0) represent the recovered/removed and infected fractions of the population at time t = 0. For small ɛ ≪ 1 one can write I(0) ≈ (1 − k)ɛ. In turn, if two initial values are known at t = 0, the basic reproduction number is given by the initial conditions. We prefer to treat the inverse basic reproduction number k and ɛ as variables, so that R(0) and I(0) are determined by this set. Any other choice of two variables from the set k, ɛ, I(0), and R(0) would work equally well, as the two remaining ones are then given by equation (12).

It it important to realize that the all-time case (I) has solutions only for a limited range of k values. While there is no recovery when k = 0 due to equation (7), for k > 1 the number of infected would drop (proven below) rather than grow from the boundary value I(−∞) = 0, which is not possible as I ∈ [0, 1]. For k = 1, I(0) = f₁(ɛ) = 1 −ɛ− e^−ɛ is negative for any ɛ > 0, and simply remains at its boundary value I(−∞) = I(0) for ɛ = 0, and thus I(t) = 0 at all times for k = 1. Finally, because I(0) ⩾ 0 must hold, not all values in k ∈ (0, 1) are compatible with the initial condition (11) at time t = 0. The requirement I(0) ⩾ 0 is equivalent with k ⩽ k_max with

This is why the remaining meaningful range of k values is k ∈ (0, k_max) for case (I). As the SIR model is assumed to hold at all times, its solution allows to calculate all fractions at times before and after t = 0. These two features make case (I) qualitatively very different from case (II).

In contrast, for the semi-time case (II) treated by Kermack and McKendrick [5] and in part B of this work, the SIR model is not assumed to hold at all times, but only at times t ⩾ 0, the time for which feed data may be available. In that case the model cannot be used to calculate fractions at times prior to t = 0, as this would under many circumstances lead to S > 1 or I ∉ [0, 1] at times prior to the 'observation' time t = 0. Because the boundary conditions of case (I) must not be respected anymore in that case, one can use arbitrary initial conditions incompatible with equation (12). The model in this setup has therefore three independent parameters such as k, S(0) and I(0), while the remaining R(0) is given by equation (1). If one chooses the initial conditions to satisfy equation (12), we are back at case (I). Case (II) therefore chooses initial conditions incompatible with equation (12) such as S(0) = 1 − and I(0) = , implying R(0) = 0. In that case one has moreover the freedom to choose k > 1, as it does not automatically lead to I ∉ [0, 1].

Case (II) is thus more flexible, but cannot be used to adequately describe the past. While case (I) can be considered as the solution to the SIR model over the whole time domain, case (II) can be used to study the future and the case k > 1 of relevance toward the end of an epidemics. In this work we are going to study case (I) and have to thus assume k ∈ (0, 1) and ɛ ⩾ 0.

2.3. Three important remarks

First, as has been noted before [9], in the three dynamical equations (2), (3) and (6) there is migration from the S-compartment (susceptible) to the I-compartment (infected) at a rate proportional to SI, and removal from the I-compartment to the R-compartment (recovered, dead or isolated) at a rate proportional to I; there is no exit from the R-compartment and no entry into the S-compartment.

Secondly, it is necessary to start initially at time t = 0 with at least one infected person, i.e., I(0) = f_k (ɛ) > 0 (see the second initial condition (12)), as the dynamical equation (2) implies for the initial change if I(t = 0) = 0. Nothing can grow out of nothing; a situation very similar to kinetic plasma instabilities where seed electromagnetic fluctuations, often from spontaneous emission [27], are needed for starting the instability. Keeping exact tract of the initial conditions during the analysis will prove to be essential to avoid mathematical singularities in the analysis.

Thirdly, the first term on the right-hand side of the dynamical equation (2) denotes the newly infected population fraction, i.e., the differential rate of the total fraction J(t) of persons that have ever been infected, i.e.

where the second equality results from equation (3). It is this rate that can be measured, is usually reported by health agencies along with the cumulative fraction , and that has been modeled by the Gaussian time evolution in our earlier work [3, 4]. Because S(0), I(0) are usually not reported as they cannot be measured directly, we will show how to replace SIR parameters k and ɛ by an initial condition for the typically available J(0) and .

Throughout this work we use the bar-notation if we approximate a quantity q.

3.1. Reduced time

We are here assuming that the infection and recovery rates have the same time dependence, so that their semipositive ratio is time-independent:

corresponding to a SIR-basic reproduction number r₀ > 1. We emphasize that this special case includes the standard case used by most analysis before that the infection and recovery rates are constants with respect to time. Equation (16) still allows us to take into account an arbitrary time-dependence of the infection rate a(t) which, of course, then is identical to the time dependence of the recovery rate μ(t) due to assumption (16). Then it is convenient to introduce for arbitrary time dependence of the infection rate a(t) the new dimensionless time variable τ with τ(0) = 0 via

so that

Using the reduced time τ, equations (2), (6), and (7) simplify to

and

respectively. In equation (19) (d ln  S/dτ)_τ=0 = −f_k (ɛ) in order to meet the initial condition (12) for I(τ = 0) = f_k (ɛ).

Equation (19) includes the well-known threshold [5, 9] value k which determines the so-called J-shape or peak-shape [9] of the pandemic wave. Obviously, equation (19) indicates that a peaked pandemic wave cannot emerge if the ratio k > 1: as S ⩽ 1 the right-hand side of equation (19) in this case is always negative, so that with positive I(τ) the time derivative dI/dτ < 0, so that the infection rate decreases with time. An epidemic wave emerges if the ratio k < 1; the case k = 1 we already discussed in section 2.2.

The boundary conditions S(t) = 1 and R(t) = 0 for t = −∞ take over to S(τ = −∞) = 1 and R(τ = −∞) = 0 for k > 0, because a positive k implies μ(t) > 0 for all t, and thus also a(t) > 0 for all t. The solution to equation (20) is thus formally equivalent to the relationship we had between R(t) and S(t),

3.2. Analytical solution

Since the original pioneering work [5] the procedure to obtain analytical solutions is similar in later work [9, 28] and also here: either (1) one expresses two of the interesting variables S(t), I(t) and R(t) in terms of the third one as done in reference [5], or (2) one expresses all three variables in terms of a suitable chosen function as done in reference [28] and here. In both cases one uses equation (1) to calculate the solution. Inserting equations (19) and (21), equation (1) reads

We emphasize that this equation fulfills the initial conditions (11) for τ = 0. For given rates a(t) and k = μ₀/a₀ the solution of equation (22) yields the temporal evolution of S(t). Equation (22) can be written as

or

which apart from a different notation corresponds to equations (23) and (27) of Harko et al [28]. The solution to the dimensional version of equation (24) had been expanded into a series by Barlow and Weinstein [15], who realized that the convergence radius is rather low and suggested to look for solutions with correct asymptotic behavior, as we will do here.

To this end we follow a slightly different but equivalent approach starting out from equation (22). We introduce the quantity G(τ) by

with the initial condition G(τ = 0) = ɛ. According to equations (19) and (21) the function G (25) determines

The dynamical equation (22) in terms of G(τ) reads

using the function f_k defined in (13). Integrating

over τ then provides (figure 1)

where the integration constant was determined from the initial condition G(τ = 0) = ɛ. The integrand vanishes at x = 0 and at certain x = G_∞ to be discussed below. The whole τ ∈ [−∞, ∞] range is thus captured by equation (29) upon varying G between zero and G_∞, and vice versa, G(τ) ∈ [0, G_∞] determines the range of G.

Zoom In Zoom Out Reset image size

The solution (29) generalizes the known analytical solutions in the literature [5, 9, 28] as it holds for arbitrary time-dependence of the infection rate a(t). The mentioned known solutions can be reproduced with equation (29) by setting τ = a₀ t on its left-hand side resulting from a constant injection rate a₀. The parametric solution presented in the original SIR work [5] assumed G ≪ 1 to simplify the analysis. We will see below that G_∞ ≪ 1 holds only within a very narrow range of k values close to k = 1. The contributions within the regime G ∈ [0, ɛ] are usually not considered in numerical schemes while they can be used to verify the proper boundary conditions stated at the very beginning. Harko et al [28] did show results for t > 0 only, as some of the fractions would get negative or exceed unity at times prior to t = 0, using their analytic solution. We will present an equivalent formulation with τ(J) of the solution (29) in equation (41).

As NPIs during the pandemic wave generate a time-changing infection rate a(t), the generalized solution (29) is highly valuable to assess quantitatively the effect of the NPIs on the time evolution of the disease.

3.3. Maximum value G_∞

The solution (29) indicates that the maximum value G_∞ = G(τ = ∞) of the function G is attained when the denominator of the integrand vanishes, i.e.

or

which is of the form (G1) with c = 1, a₀ = −k and r = 1/k. Consequently, according to equation (G2) we obtain (figure 2)

where W₀ is the principal solution of Lambert's equation z = W(z)e^W(z) (discussed in appendix G) and α is given in terms of k as

The argument α of the Lambert's W₀ function in equation (32) is negative for all k (figure 15(b) in appendix G). Accordingly, W₀ ∈ [−1, 0] is negative as well (figure 15(a)), while G_∞ ⩾ 0. In figure 2 we show the maximum G_∞ as a function of k. It can be seen that G_∞ has values significantly greater than unity for k ≪ 1, while it vanishes at k = 1. The value G_∞ = 1 is attained at k = 1 − 1/e ≈ 0.63. The dashed black line shows the asymptotic behavior of G_∞ ∼ 1/k at small k ≪ 1. Significant deviations from the asymptotic behavior set in at k ≈ 1/3. According to these considerations, G_∞ is well approximated by (figure 2)

The relative error of this approximant is below 2.2% for k ⩽ 0.63; its absolute error is below 0.014 for k ⩾ 0.63, while it is exact at k = 1 and exhibits the correct asymptotic behavior for k → 0.

Zoom In Zoom Out Reset image size

Similarly, the non-principal solution of Lambert's equation sets the lower bound G_−∞ = G(τ = −∞) with G_−∞ = k⁻¹ + W₋₁(α) = 0. The lower bound is identical zero because Lambert's equations is trivially solved by W₋₁ = −1/k for α given by (33). To summarize, G ∈ [0, G_∞], compatible with τ ∈ [−∞, ∞] and S ∈ [0, 1]. This compatibility is established as a result of the proper boundary conditions.

3.4. Differential and cumulative rates of newly infected persons

From the invariant dt = j(τ)dτ we obtain with equations (25), (27), and equation (17) in the form dτ/dt = a(t) for differential rate of newly infected persons (figure 3)

with the initial value j(τ = 0) = f_k (ɛ)e^−ɛ . The corresponding cumulative distribution picks up all the newly infected individuals, not only the ones that are going to occur at t > 0, i.e., one has (figure 4)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the far past, both quantities j and J emanated from zero, because G_−∞ = 0. As indicated we may regard these j and J as functions of G(τ).

Equation (19) in the form dS/dτ = −IS, and dI/dτ = SI − kI imply with j = SI from equation (35)

A maximum of  ln  j and j thus occurs when S − I = k, while (25) and (28) can be used to write S − I in terms of G

Hence equation (38) can be written as

Having introduced j and J = 1 − e^−G , it is worth mentioning that the solution (29) can also be expressed in terms of J as follows. Since dG/dJ = (1 − J)⁻¹,

with J(ɛ) = 1 − e^−ɛ . The range of J values follows from the range of G values, and is thus given by , in light of equation (31). With τ(J) at hand, the remaining quantities of the SIR model are obtained immediately via S = 1 − J, R = −k ln(1 − J), I = J − R, and j = SI,

3.5. SIR parameters

The J thus captures all infections up to τ, including those that occurred prior τ = 0, and the usually reported J(0) = 1 − S(0) = 1 − e^−ɛ provides us with S(0) and ɛ, i.e.,

This possibility to determine a parameter does not exist for case (II) where such quantity J is not available from the model equations and can only be manually entered as a new parameter. Only for case (I) it is directly related to S(τ), as we have just shown. Having determined ɛ from J(0), and assuming τ = a₀ t with a constant a₀ = a(t = 0), there are two remaining parameters a₀ and k of the SIR model, that we now obtain from the measurable J(0) and its derivatives and at the initial t = 0. To this end we write down expressions for and and solve them for k and a₀. The first relationship is obtained from given by equation (36), evaluated at t = 0. This yields with f_k (ɛ) = 1 − kɛ − e^−ɛ that

The second relationship follows from , using equations (27) and (36). As dj/dτ = (dj/dG)f_k (G) we can proceed and express in terms of G(0) = ɛ, k, and a₀,

where the 2nd line can be used if ɛ ≪ 1. Under such conditions

Upon replacing k from equation (43) in equation (45) one can write down an explicit, but lengthy, expression for a₀ in terms of J(0), and . We have thus shown how to obtain all three SIR parameters from the cumulative fraction of newly infected persons J and its derivatives at t = 0. In other words, the parameters are obtained from a quadratic fit to the reported J(t), taking into account only data in the vicinity of t = 0. Depending on the available data, other expressions presented in this work may be preferred to extract the coefficients, such as equation (68).

3.6. Final values

The knowledge of G_∞ from equation (32) and the solution (29) allow us already to derive a number of important results. We first consider the final values at τ = ∞ of (figure 2)

where we have used equations (19), (21) and (25). As a consistency check we note that equation (1) is fulfilled here, i.e. S_∞ + R_∞ + I_∞ = 1. Inserting G_∞ from equation (32) we obtain with the help of equation (32)

which both are determined solely by the inverse basic reproduction number k = 1/r₀. In figure 2 we plot the fractions R_∞ and S_∞ from equation (48) as a function of r₀ > 1 along with their approximants that follow immediately follow from equations (34) and (46).

If the final values were known, one can calculate the corresponding r₀ upon inverting the relationship between S_∞ and r₀. The result is r₀ =  ln(S_∞)/(S_∞− 1). The canonical values R_∞ = 2/3 and S_∞ = 1/3 are thus reached for r₀ = 3  ln(3)/2 ≈ 1.65. The final susceptible fraction S_∞ decreases with increasing r₀, starting with S_∞ = 1 at r₀ = 1, and reaching S_∞ = 0 as r₀ grows.

For these values of the differential and cumulative rates of newly infected persons after infinite time we obtain with equations (35) and (46)

J_∞ thus coincides with R_∞ (figure 2). Its approximant based on equation (34) is

With increasing values of k the cumulative fraction of infected persons decreases from 1 at k = 0 to 0 at k = 1.

3.7. Bell-shaped differential rate. GM-like solution of the SIR

Equation (36) can also be used, even without the explicit inversion of equation (29), to derive generally valid expressions for the time τ₀ of maximum, the maximum level j_max and the dimensionless width ω of the differential rate j, which correspond to the three important parameters characterizing a Gaussian differential rate,

although here the width ω may not be a τ-independent constant. Assuming a constant ω, equation (51) is known as the GM for the time evolution of the daily number of new infections (or also deaths) [4].

3.7.1. Exact peak amplitude

The maximum of the differential rate (36) occurs when the derivative dj/dτ in equation (40) vanishes providing

directly from equations (38) and (39). Writing this as

makes the equation of the form (G1) with c = 1, a₀ = −k/2 and r = (1 + k)/k. The analytic solution of equation (52) is the non-principal solution W₋₁ of Lambert's equation,

with α from (33). If W in (54) were the principal solution W₀ of Lambert's equation, with the property W₀(α) ∈ [−1, 0] and W₀(α₀) ∈ [−0.406, 0] for k ∈ [0, 1] (appendix G), then G₀ = G_∞ + 1 + W₀(α₀) − W₀(α) would exceed G_∞. The W in (54) is therefore the non-principal solution W₋₁ of Lambert's equation with the property W₋₁(α₀) ⩽ − 2 for all k ∈ [0, 1], as shown in appendix G. A simple approximant for G₀ is

It is exact at k = 0 and k = 1 and has a maximum absolute error of 0.008.

Inserting equation (52) in equation (36) yields an analytic expression for the maximum value

where the 1st, 2nd, and 3rd line are obtained with the help of (36), (54), (53), respectively. Since W₋₁(α₀) ⩽ − 2 for all k ∈ [0, 1], the j_max is positive. For k = 1 the W₋₁(α₀) = −2 and j_max vanishes. For k → 0 the W₋₁(α₀) diverges to −∞ (appendix G), but the k² in front causes j_max to approach 1/4. Note that j_max (58) is not affected by ɛ but is solely determined by the inverse basic reproduction number k. This is reflected by the results shown in figure 3: regardless of ɛ the maximum value j_max remains unchanged for given k. A simple approximant for j_max is provided by equation (77).

We recall that the maximum j_max in the measurable daily number of new infections is completely different from the maximum population fraction I_max of the infected compartment. Because I(τ) = f_k (G) it achieves its maximum at I_max = 1 − k + k  ln  k according to equation (A2).

3.7.2. Peak time

The time, at which j_max is reached, need not be positive. The maximum may have occurred already, depending on the initial condition and the value for the inverse basic reproduction number k. According to solution (29) the value of G₀(k) corresponds to the time of maximum τ₀(k) where G₀ and τ are given by equations (54) and (29). An explicit expression for τ₀ in terms of k and ɛ is given upon inserting G₀ from equation (54) into equation (86). This yields (figure 5)

with the very useful abbreviations

and (figure 6(a))

where 1 + W₀(α) is identical to the crossover x* from equation (A11). We can use this κ with the feature κ ∈ (0, 1) to rewrite the exact expression (32) as

Figure 3 indicates that the peak time τ₀ varies inversely with ɛ, confirming equation (59).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.7.3. Exact behavior in the vicinity of the maximum

With the help of equations (28) and (40) one has

so that a 2nd order Taylor expansion of  ln  j(τ) around its maximum

can be evaluated, where ln  j_max is given by equation (58). Making use of as given by (53), equation (63) evaluated at τ = τ₀, or equally, G = G₀, simplifies to

with (figure 7(b))

We thus have for the Taylor expansion of  ln  j(τ) around its maximum value at τ₀ to second order, with

with ω₀ given by equation (66).

Zoom In Zoom Out Reset image size

3.7.4. Exact asymptotic behavior of the SIR model

Next we show that the rate j(τ) exhibits exponential behavior at early and late times with respect to the peak time and evaluate the growth coefficients.

Making use of the limiting values G(τ = −∞) = 0 and G(τ = ∞) = G_∞ in equation (40) we find the growth coefficients characterizing exponential growth and decrease at early and late times,

with κ from equation (61) and τ₀ from equation (59). Note that the exponent in (68) is positive, and the exponent in (69) is negative since κ > 0 for all k ∈ (0, 1).

To calculate the important prefactors in these expressions we start from equation (B2). At small times, where g = G/G_∞ approaches zero, we use the Taylor expansion  ln(1 − g) = −g to obtain

for g ≪ 1 as the solution of equation (B2). Since the argument of Lambert's principal solution W₀ gets small in this limit, cf equation (G10), because j = (1 − k)G + O(G²), and with the help of τ₀ (59) and A(x) (60) we obtain for the regime of early times τ ≪ τ₀

In the opposite limit, where g = G/G_∞ approaches unity, we can use the Taylor expansion  ln(g) = g − 1 + O(1 − g)² to obtain

for 1 − g ≪ 1 as the solution of equation (B2). Since the argument of Lambert's function W gets again small in this limit, because , we obtain for the regime at late times τ ≫ τ₀

We have thus shown that the initial condition S(0) enters only the peak time τ₀ and that the prefactors j_early and j_late characterizing the amplitudes of the asymptotic behaviors depend on k alone. Their dependence on k is shown in figure 8. While j_early reaches its maximum at k → 0, where the SIR reproduces the SI model (section 3.6), j_late becomes very small in this limit. As a matter of fact, the asymptotic behavior of the SI model, corresponding to a value unity in figure 8, is all that can be observed up to huge times τ ∼ k⁻¹ (denoted as τ_c in appendix F), as not only the prefactor j_late becomes extremely small in the limit k → 0, but also the contribution j_late/κ(1 − k) from the asymptotic exponential to the cumulative J.

Zoom In Zoom Out Reset image size

These expressions (68), (69), (71) and (73) exactly capture the short and long time behavior of j(τ) shown in figure 3. The coefficients ɛ and k can thus also be read off from the regime of exponential growth, if a₀ is already known. It is worthwhile noticing that the initial condition S(0) does only affect the peak time, but not the prefactors j_early and j_late in the asymptotic behavior of j (figure 9).

Zoom In Zoom Out Reset image size

3.7.5. Gaussian width

While j is not a perfect Gaussian, as it is not perfectly symmetric (to be discussed below), we can calculate a width via three routes. Route (A) fits j(τ) by the Gaussian function (51) with constant w. Route (B) determines w by the known GM relationship using our above results for J_∞ and j_max. The third route is to estimate a width from the behavior in the vicinity of the peak times. The latter leads to width ω₀ (66). The former two routes give rise to very compatible widths ω as a function of k (all widths shown in figure 7(b)).

Making use of our above results (49), (46) and (32), we can thus express a τ-independent characteristic width of j(τ) analytically as

where α₀ = 2α/e and α expressed in terms of k in (33), and where we have used G₀ from equation (54). Note, that the width is unaffected by ɛ as well, it is solely determined by the inverse basic reproduction number k. Taking into account all limiting behaviors our approximant for ω is given by

implying, with the help of our approximant (34)

Both are compared with the exact results in figure 7. The two frequencies ω₀ and ω are actually very comparable which lead us to the GM-like approximant to be discussed next. We recall that the dimensional for a constant a(t), and that the dimensional width is then w = ω/a₀.

3.7.6. GM-like approximant for j(τ)

We have by now calculated the behavior of  ln  j(τ) around its maximum including peak height (denoted as j_T (τ) in equation (67)), and the asymptotic behavior of j(τ) at small and large τ. Using a continuity requirement we can construct an approximant for j(τ) that captures the behavior of j(τ) qualitatively at all times. The continuity requirement is, in light of equations (68) and (69),

where j_early and j_late are given by equations (71), (73), and where the two times τ_early and τ_late have to be determined by the continuity conditions. The condition (78), with j_T (τ) from equation (67), yields

Similarly, condition (79) yields

The full GM-like approximant for j(τ) therefore reads (figure 10)

with j_early and j_late given by equations (71) and (73), κ defined by equation (61), τ_late and τ_late according to equations (80) and (82), ω₀ specified by equation (66), and τ₀ denoting the peak time τ₀ = τ(G₀) from equation (59). The classical GM [4] is a special case of the GM-like approximant for τ_early = −∞ and τ_late → ∞.

Zoom In Zoom Out Reset image size

Not under all conditions does the central Gaussian part meet the asymptotic exponential tail. Both D_early and D_late can become complex-valued as soon as the argument of the square root becomes negative. This is the case for D_late if τ₀ is sufficiently large and negative. Under such conditions the relevant crossover times are the real parts of τ_early and τ_late, and the asymptotic branches are only shown up to the level that is reached by j_T (τ).

Our GM-like approximant (84) suggests that the GM is a suitable approximation for times between τ_early and τ_late. The ω from equation (74) and ω₀ from equation (66) are found to be very similar, which implies that J_∞ calculated for the GM can serve as an estimate for the exact J_∞, as done in reference [4]. The difference between the GM and the SIR can also be visualized by eliminating time, and plotting j/j_max versus J/J_∞ (figure 11).

Zoom In Zoom Out Reset image size

The value of τ_early can be most easily determined from the differential and cumulative doubling times [29] of the GML-approximant (84), as both are constants below t_early. t_early coincides with the time when the monitored doubling times start to increase with time.

Similarly, τ_late can be determined when the differential decay time (although this is seldomly reported) or the cumulative doubling time become constants at late times of a pandemic wave.

The cumulative

as well as the final corresponding to the GM-like then follow by integration. Analytic expressions are provided in appendix H.

In order to derive explicit analytical expressions G(τ) it is necessary to calculate, at least approximately, the remaining integral on the right-hand side of equation (29). Whereas in the references [9, 28] the analytical evaluation of this remaining integral is not detailed, we recall that in the pioneering work [5] this was done for small values of G less than unity by expanding the exponential function in the denominator to second order in x about x = 0, i.e. e^−x ≃ 1 − x + x²/2. However, this approximation cannot be applied here for all values of G as figure 2 indicates that G_∞ attains values significantly greater than unity. Therefore more generally, without referring to the regime x ≪ 1, an analytic approximant for τ can be derived, as shown in appendix A

that is valid for all k ∈ (0, 1) and all ɛ ≪ 1. The κ is defined in terms of k by equation (61). We note that equation (86) can alternatively be written as with A defined by equation (60). Obviously, is fulfilled for this approximant. Upon replacing G by G₀ in (86), we arrive at the expression (59) for the peak time. In figure 1 we compared the numerically evaluated integral (29) with the approximation (86) for various k. We notice the excellent agreement.

Even more important in practice, we were able to invert the relationship (86) to obtain an approximant . This is done in appendix (B). There we use this approximant to confirm that it exhibits the correct asymptotic behaviors, that it reproduces the exact j_max , and shares the peak time τ₀ with exactly. To be specific we obtain

where the k-dependent parameters b_1,2 and c_1,2 (figure 6(b)) are specified in appendix B, τ₀ given by equation (59), and W₀ is the principal solution of Lambert's equation (appendix G). With at hand, one can express all results in terms of τ, and thus also in terms of time t—for arbitrary a(t)—without preparing a parametric plot. The approximant (87) is implemented in appendix J.

As we have demonstrated, J_∞ and all other final values are insensitive to the precise time-dependency of a(t). It is possible to flatten the curve of daily infections by interventions that reduce a(t_max) at peak time t_max given by τ₀ = τ(t_max). Here is the maximum population fraction that is newly infected within a single day, if we choose day as the time unit in which we also express the infection rate a(t).

While J(t) = J(τ) is invariant, the differential rates are related as . Inserting τ(t) for an arbitrarily chosen a(t) into the expression for j = f_k (G)e^−G with from equation (87), the daily number of new infections is available analytically as a function of time.

Starting from equations (36) and (37) one has G = − ln(1 − J) and S = 1 − J, so that we can write down a relationship between the differential j(τ) and the cumulative J(τ) number of infected persons,

The same relationship follows from the derivative of equation (41) with respect to τ. Since , the infection rate a(t) can be extracted from the measured and J(t) via the right-hand side of equation (88), divided by , as long as k is time-independent.

We demonstrated that the exact parametric solution S(t), I(t) and R(t) to the SIR model with time-dependent a(t), time-independent k and arbitrary initial conditions is determined by the function G(τ) uniquely specified equation (29). The t-dependency of G is given by τ(t) from equation (17). With G_t = G(τ(t)) as a function of t therefore at hand,

and I(t) = 1 − S(t) − R(t) solve the SIR equations at all times, and satisfy the boundary and initial conditions exactly. The measurable differential rate of newly infected and its cumulative counterpart J are then available from and J(G_t ) = 1 − S(t), according to equations (36) and (37). Alternatively, one can go over the J(τ) route instead of G(τ) as explained in section 3.4.

While the SIR model is usually solved numerically, and the relationship between G and τ given by equation (29) can also be calculated numerically, we provide analytic approximants (87) as well as for the inverse (86) that both have the correct asymptotic behaviors and may be use instead of the exact result, as they are sufficiently accurate for any practical purpose. Using this approximant (87), that works for any k and any initial condition, the usually reported daily number of new infections we have expressed explicitly as a function of time t for an arbitrarily time-dependent infection rate a(t) in section 5. The presented solution of the all-time SIR model can thus be used to replace the numerical solution of the SIR model. Having an analytical solution should be also advantageous to understand the effect of interventions reflected by a(t) much easier, and in a more transparent fashion.

The whole shape of the differential rate j depends on k only. The shape is characterized by the maximum j_max (58), width ω (74), asymptotic exponents 1 − k and − κ(1 − k) as well as their prefactors j_early (71) and j_late (73), while the initial condition affects the peak time τ₀ (59). All these quantities have been determined analytically in this work for the all-time SIR model. In addition, we have provided approximants for all quantities characterizing the exact solution of the SIR model, to appreciate their qualitative behavior, or in case the Lambert functions are unavailable. For the special case of time-independent a(t) = a₀, and thus τ = a₀ t, the SIR parameters k, ɛ, and a₀ can all be determined from initial conditions J(0), , and by equations (42), (43) and (45). This is important because the I and R fractions are usually not reported at a certain time to provide an initial condition.

Our derivation not only generalizes existing solutions, but may appear very compact compared with existing ones, while it includes them as a special case. We have relaxed the assumptions of t ⩾ 0 and G ≪ 1 from previous approaches, and present the unique solution that captures consistently both the future and the past, if the SIR model is assumed to hold in the past as well. The classical assumption G ≪ 1 we have shown to fail except at inverse basic reproduction numbers close to unity.

The evolution of the usually reported differential rates assumes a bell-shaped curve characterized by a peak time τ₀, height j_max, and width ω. We have provided analytic expressions for all these parameters by equations (59), (58), and (74), respectively. They involve the two solutions W₀ and W₋₁ of Lambert's equation, usually available in scientific software. To appreciate the validity of our arguments and derivations, we have collected the required exact features of W₀ and W₋₁ in appendix G.

We have inspected the limiting case of k → 0 of the SIR model to find that the SI model, corresponding to k = 0, does not trivially emerge from the SIR model in this limit. While the asymptotic behavior of the SIR model is qualitatively different from the SI model, the quantitative difference between them still decreases with decreasing k. Insofar is the SI model compatible with the limiting case of an SIR model, but at the same time asymptotically different.

The semi-time SIR model, which can be considered as the SIR model subject to initial conditions incompatible with the all-time SIR model, and thus useful to predict only the future time-evolution, will be elaborated in part B of this work.

Here we address the formal solution (29)

with f_k defined in (13), over the range of all possible values for G ∈ [0, G_∞] and k ∈ [0, 1], where G_∞ is given by (32). For the special case of G = G₀, we thus also address the peak time τ₀ = τ(G₀) introduced in section 3.7. The f_k in the denominator of the integral is semipositive, has its maximum at x_m = − ln(k), and reaches zero at G_∞. A 2nd root of f_k is located at x = 0. Accordingly, f_k has the following properties: f_k (0) = f_k (G_∞) = 0, f_k '(0) = 1 − k, f_k '(x_m ) = 0, and

For x ≪ 1 the function f_k (x) can be expanded about x = 0 as with

Likewise, for large x we expand f_k about G_∞ to obtain with

with κ from equation (61). Using these asymptotes we obtain the following approximant

Examples for the cases of k = 0.3 and k = 0.6 are shown in figure 12. Both partial integrals can be performed analytically as

Zoom In Zoom Out Reset image size

and . While τ₁ dominates the integral at small G, at G = ɛ and beyond, τ₂ is responsible for the sharp increase at large G, when G approaches G_∞ (see figure 1). With A(x) defined by equation (60), can alternatively be written as

The exact asymptotic behaviors at small and large G are

where the intermediate expressions are also required, e.g., to calculate j_late in equation (73).

As already noted, h₁(x) and h₂(x) dominate the behavior of the integrand in distinct x regimes. To specify the regimes, we solve h₁(x*) = h₂(x*) to obtain

with the help of (32), and where we recall that W₀(α) ∈ [−1, 0) is negative for all k ∈ [0, 1]. For x < x* the dominating contribution is h₁(x). While x* = 1 for k = 0 and x* ≈ 1 for k < 0.2, it depends almost linearly on k for larger k and reaches x* = 0 for k = 1. From equation (A6) we infer

and G₀ given by equation (54) we can derive from equation (A11) the inequality

for all k, while G₀ = x* at k = 1 follows from the mentioned properties of Lambert's functions.

Starting from the high quality approximant derived in the previous section, we are looking for an explicit expression for the inverse function, . Unfortunately, cannot be inverted analytically over the whole G domain, but we are going to find approximants for G separately over the two disjoint intervals and so that

with g₁, g₂ given in terms of k and ɛ by equations (B4), (B6) to be derived in this section, and compared with the exact G(τ) in figure 13. We recall that G_∞ is known in terms of k by equation (32). The two partial approximants meet exactly at τ = τ₀ specified by equation (59).

Zoom In Zoom Out Reset image size

To derive all quantities appearing in equation (B1) we start from in terms of G, provided by equations (A6) and (A7). Upon introducing g = G/G_∞ and the known g₀ = G₀/G_∞, using the abbreviation κ already introduced by equation (61), the governing equation for g is written

involving the g-independent and semipositive constant A(ɛ) defined by equation (60). Equation (B2) leads to the exact transcendental equation

that cannot be solved analytically. To proceed we approximate, for the two cases of g > g₀ and g < g₀ separately, corresponding to τ ⩽ τ₀ and τ ⩾ τ₀, respectively, the smaller of the two logarithmic terms. In each case we make sure that the approximation is correct in all limiting cases, so that not only the asymptotic behavior of G and j will be correctly captured, but also g₁ and g₂ will coincide at τ = τ₀.

B.1. Approximant

To obtain g₁ we approximate  ln(1 − g) ≃ −C₁ g in equation (B2) with a constant C₁ in such a way that the approximation is exact at both g = 0 and g = g₀. This yields . If we furthermore introduce x via g = e^−x , the corresponding equation for x is of the form (G1) with c = 1, r = A(ɛ)/κ − (1 − k)τ, and a₀ = κ/C₁. Its solution is thus x = r + W₀(e^−r /a₀), involving W₀ and not W₋₁ because the argument e^−r /a₀ ⩾ 0 (figure 15). With x at hand, we have found g₁ =   ln  x which can be written with the help of Lambert's equation in two equivalent ways

both for a different purpose. In equation (B4) we have introduced the abbreviations b₁ = κ/C₁ = −κg₀/ ln(1 − g₀),

with g₀ = G₀/G_∞, so that x₁ carries the dependency on time, and is unity at τ = τ₀ with τ₀ according to (59), G₀, κ and A explicitly given by equations (54), (61), and (60). This expression for τ₀ is just another way of writing the expression (59) for τ₀ we obtained from in the previous section. The equivalence follows already from equation (B2). The right-hand side of equation (B4) proves that g₁(τ₀) = g₀.

B.2. Approximant

For the opposite case of large g close to unity. We use the same strategy as before and now approximate  ln(g) ≃ C₂(g − 1) with a constant C₂ in such a way that the approximation is exact at both g = 1 and g = g₀. This yields . If we proceed in a way analogous to the case of g₁, we obtain

Here we have introduced the abbreviations ,

so that x₂ carries the dependency on time, is unity at τ = τ₀, where τ₀ was already given by equation (59).

The two half approximants (B4) and (B6) (figure 13) perfectly match at G₀ by construction, i.e., by the proper choice of the approximants for  ln(1 − y) and  ln  y. The calculations leading to the present form of the approximants had to make use of properties of Lambert's function, in particular the identity W(z ln  z) =   ln  z from equation (G6) for any positive z. Note that the dependency of on ɛ is adsorbed by x₁ and x₂, while all other coefficients are positive and depend on k only. The coefficients are not independent of each other, as

holds. To summarize, in this section we have provided an explicit expression for the approximant in terms of k and ɛ, and also an approximant for the peak time τ₀ in terms of the same parameters. Because the assumptions were asymptotically exact, the approximant has the features

as expected, and is compared with the exact G(τ) in figure 13. To arrive at equation (B9), we used property (G10) of the Lambert's function. Approximants and are compared with each other in figure 14.

Zoom In Zoom Out Reset image size

S(τ) had been defined by equation (25). Using the identity W(z) =  ln(z) −  ln  W(z) =   ln[z/W(z)] implies aW(z) =  ln[z/W(z)]^a and e^−aW(z) = [W(z)/z]^a . We thus have, starting from the result of the previous section,

where captures the regime before, and the regime after the peak time τ₀. As these approximants have correct asymptotic behavior, the exact asymptotic behavior of S(τ) follows from these expressions with equation (G10).

Since R = kG according to equation (26) and S + I + R = 1, we can write immediately down our approximants for R(τ) and I(τ)

with and given by equations (B1) and (C2).

The differential rate j(τ) had been defined by equation (35) or the equivalent equation (36). With from equation (B1) at hand we can calculate the approximant for via equation (36).

with τ₀ according to equation (59). with the approximants for and in terms of already written down in equations (D3), (C1), and (C2). Obviously, j₁ meets j₂ at τ = τ₀, because G₁(τ₀) = G₂(τ₀). Note that this is not the case if equation (35) is used, as the slopes of and are not exactly identical at τ = τ₀. The prefactors of this approximant in the asymptotic limits of early and late times differ from the exact equations (71) and (73) only by a factor C₁ and C₂, respectively, introduced in section B. Both factors approach unity for small ɛ, and start to deviate from unity only at large ɛ ≳ 0.2.

The maximum of the approximant coincides with the exact j_max , because we have determined using equation (36), and because the definition of G₀ implies , which remains unchanged, as the maximum of the approximant occurs at τ = τ₀, where G = G₀.

E.1. Remark: non-monotonous

If we calculate the approximant starting from equation (35), using the notation from appendix B, with dx₁/dτ = (1 − k)x₁, and the property (G8) of Lambert's function, this gives

Similarly, with dx₂/dτ = −κ(1 − k)x₂,

At τ = τ₀, these expressions simplify with x₁ = x₂ = 1, g₁(τ₀) = g₂(τ₀) = g₀, . We thus have

If we now make use of relationships from appendix B, i.e., replace W₀(c₁) = −κ⁻¹  ln(1 − g₀), W₀(c₂) = −κ  ln  g₀, and κ by its definition (61), we see that this ratio is unity for a single k, but not all k ∈ (0, 1). This is because is generally non-monotonic.

For the extremal case of k = 0 (so-called SI model), which is shown here to emerge nontrivially from the SIR model in the limit k → 0, the differential equation (27) for G(τ) with G(0) = ɛ is solved by

hence and

since R_SI = 0 at all times for k = 0. The differential rate is given by G_SI(τ) via equation (36) or via equation (35), or from dJ_SI/dτ, as

The same result is obtained with j_SI(τ) = I_SI(τ)S_SI(τ) using the solution of the logistic equation, mentioned in section 2.2. Its maximum with amplitude thus occurs at , the asymptotic behavior is given by and , and the GM-like behavior in the vicinity of follows from the Taylor-expansion as , corresponding to ω₀ = 2. Within the SI model, every person gets infected in the course of time, , but nobody will ever get uninfected again (through recovery/removal).

Importantly, while some quantities like , , obviously agree with the values one obtains for the SIR in the limit k → 0, the exponential decay of j_SI (F3) of the logistic SI model (k = 0) is qualitatively different from the SIR model at times after peak time, as the SIR exponent −κ(1 − k) tends to reach zero for k → 0. In fact, the deviations between SI and the k → 0 limit of the SIR get smaller with decreasing k. This can be seen by inspecting τ in the regime where it is large, and G close to G_∞. For small k one has κ ≃ k and κ(1 − k) ≈ k, as well as G_∞≃ k⁻¹ and thus, according to equation (A7),

as k approaches zero. This relationship can be inverted to write G ≃ k⁻¹(1 − e^−kτ ). Since e^−kτ = 1 − kτ + O(k²) for any finite τ, the G approaches τ, and j thus approaches, albeit very slowly with decreasing k, the limiting , in agreement with the SI model (F3). To be more precise, the SIR approaches the SI in the sense

while for any small but finite k, the SIR j still enters the regime j ∝ e^{−κ(1−k)τ} ≃ e^−kτ at some finite crossover τ = τ_c < ∞ that can be estimated by equating this asymptotics with equation (F3). The asymptotics of the SIR at large times is thus different from the SI, but the numerical difference between the SI and SIR results still decreases with decreasing k. To estimate τ_c in the limit of small k, one can safely use the mentioned approximations G_∞ ≃ k⁻¹ and κ(1 − k) ≃ k. By equating the two differential rates (69) with (73) and (F3) the governing equation for the crossover τ_c is thus , or equivalently,

where we have used another relationship, k(1 − k) ≃ k, also valid for small k. Further assuming k ≪ɛ≪ 1, this latter expression simplifies to

The crossover τ_c thus increases not only strongly with decreasing k but moreover slowly with decreasing ɛ, and this makes equation (F5) to hold.

Another signature of qualitatively different setting between the SI and SIR at k → 0 is the relationship J_∞ = R_∞ for the SIR, which breaks down at k = 0, as R does not depend on time for k = 0, and thus and , while the limiting behavior of the SIR model is  lim_k→0 R =  lim_k→0 kG_∞ = 1.

Irrespective these observations, for an arbitrary infection rate a(t), the daily number of new infections as a function of time t is given for the SI model by equation (F3) with .

Lambert's W function [30, 31] with several applications in physics [32–37] solves the equation z = W(z)e^W(z). It can be used to calculate the solution of transcendental equations of the form

in explicit form as

The two real-valued solutions are denoted by W₀ (principal) and W₋₁ (non-principal) as shown by figure 15(a). These two functions are available in scientific software such as lambertw (Matlab^R) or ProductLog (Mathematica^R).

Zoom In Zoom Out Reset image size

For the purpose of this manuscript we are interested in W(α) and W(α₀) as a function of k ∈ [0, 1], where α = −e^−1/k /k (32) and α₀ = 2α/e (54). The two α's and their ranges are depicted as a function of k in figure 15(b). While α ∈ [−e⁻¹, 0] for any k, the α₀ ∈ [−2e⁻², 0]. The lower bound is in each case reached at k = 1. One has W₀(0) = 0, W₋₁(0) = −∞, W(−e⁻¹) = −1, W₀(−2e⁻²) ≈ −0.406, and W₋₁(−2e⁻²) = −2 as one can verify by pluggin in these values into Lambert's equation. Accordingly, for the case of α the solutions of Lambert's equation reside in the disjunct ranges W₀(α) ∈ [−1, 0] and W₋₁(α) ∈ [−∞, −1]. For α₀ the solutions reside in the disjunct ranges W₀(α₀) ∈ [−0.406, 0] and W₋₁(α₀) ∈ [−∞, −2]. These ranges are used to select the relevant solutions of Lambert's equation in (32) and (54). A special, trivial but important case is z = −e^−1/k /k = α(ɛ = 0) because this z obviously solves Lambert's equation exactly with the simple W₋₁(z) = −k⁻¹, i.e., one has

for all k. Using the same α's, Lambert's W functions have the following property

that implies G₀ ⩽ G_∞.

Proof. As is obvious from figure 15(b), the difference between α and α₀ is largest at k = 1 and ɛ = 0. For these special values W₀(α) = W₀(−e⁻¹) = −1 and W₋₁(α₀) = W₋₁(2e⁻²) = −2. The difference then changes monotonously and becomes huge in the limit k → 0 (and α → 0), as is obvious from figure 15(b). This completes the proof.

Another important inequality is

as it implies G₀ ⩾ G_−∞ = 0.

Proof. The left-hand side increases with increasing k, as α − α₀ increases as well (figure 15(b)). At k = 1, W₋₁(α) = W₋₁(−e⁻¹) = −1 and W₋₁(α₀) = W₋₁(−2e⁻²) = −2. This completes the proof.

There is another identity that we need to derive the final expression for our approximant .

holds for any z > 0. As a proof one can start from this expression for W, to obtain We^W = ( ln  z)e^ln z = z ln  z, which is Lambert's equation. A function W with the property (G6) fulfills Lambert's equation at the same time, as long as W is real-valued; equation (G6) is thus equivalent to Lambert's equation over the domain of positive arguments.

Taking the derivative of Lambert's equation z = W e^W with respect to z, Lambert's functions W₀ and W₋₁ both obey the differential equation

so that

Relationships (G8) and (G9) are required if one attempts calculating the approximant from .

For small z, Lambert's function W₀ has the expansion

implying

The cumulative corresponding to given by the GM-like equation (84) follows by integration. One has

with

Because our analytical approximants seem to reproduce the exact result so very well that it is hard to see any difference between the curves by eye, figure 16 provides some plots that show the absolute relative difference between the two.

Zoom In Zoom Out Reset image size

For the convenience of a reader, this section contains a complete Mathematica^TM notebook that creates time-dependent SIR results using our for an arbitrary infection rate a[t], specified in the first line, arbitrary k ∈ (0, k_max) and S(0) ∈ (0, 1). All results such as S(t) (code SIRS[k,S0][t]) and (code SIRj(k,S0][t]) are available analytically, and can be plotted as shown.