The theory of the SIR model is closely related to two special functions – the Lambert W function and the Wright’s Ω function, which is less well-known and can be expressed as a composition of the W function.

The Lambert W function is defined implicitly by the functional equation W z e W z ≕ z , z ∈ C (4) and as such is the simplest example of a root of an exponential polynomial. Properties of the function are surveyed in [14]. It is particularly useful in solving equations where the variable appears both in the base and the exponent, and it finds applications in various fields, such as combinatorics, physics, biochemistry, and complex analysis. Applications in physics include the Wien’s displacement law; in biochemistry – the Michaelis-Menten kinetic equation, etc. In general, the W function is multivalued. Analytically continued on the complex plane, the function has a countably infinite number of complex-valued branches. The real-valued branches are two. The principal branch is conventionally denoted by W 0 , while the non-principal real-valued branch is denoted by W − 1 . In the present work, I will slightly depart from the convention and will denote W − 1 ≡ W − and W 0 ≡ W + since we will be concerned only with the real-valued branches. Furthermore, if no subscript is written the symbol W will denote all branches of the function.

The Lambert’s function is a simple example of a function which is non-Liouvillian – it can not be expressed by a finite composition of elementary functions or their integrals [15].

Useful identities are Equations 5–9: e n W z = z W z n (5) log W + z = log ⁡ z − W + z , z > 0 (6) W + z = log z W + z , z > 1 / e (7) W − z = log z W − z , z ∈ − e − 1 , 0 (8) W n z n W z n − 1 = n W z , n > 0 , z > 0 (9)

Furthermore, the Lambert function obeys the differential equation W ′ x = e − W x 1 + W x = W x x 1 + W x

The point x = − 1 / e is a branching point of the function where the function can be presented by a convergent Puiseux series. Around the branch point, the series can be given as W p = − 1 + p − p 2 3 + 11 27 p 3 − 43 540 p 4 + ⋯ where p = ± 2 ( e x + 1 ) for W + or W − , respectively. The series converges for | p | < 2 . This relationship can be used to evaluate the function around the branch points where the usual Newton’s method converges very slowly.

The Wright’s Omega function is closely related to the W function [16]. Corless and Jeffrey [16] define the function as Ω z ≔ W K z e z , K z ≔ I m z − π 2 π (10)

The Ω function obeys the rational autonomous differential equation Ω ′ = Ω 1 + Ω (11)

This equation arises in biochemistry as the Michaelis-Menten enzyme kinetic equation. The present work uses different convention Ω ± z ≔ W ± − e z (12)

The following results are stated and proofs are repeated for convenience.

Proposition 1

The SIR system is reducible to the non-linear differential equation for the i-variable: i ′ = − g i 2 − i i ′ + i ′ 2 i (22) or the system I t ′ = − g I 2 − I I t + I t 2 I (23) I ′ = I t (24)

Proof. from the conservation law s ′ = − i ′ − r ′ = − i ′ − i g . Then from Equation 14 it holds that s = g + i ′ / i . Differentiating Equation 14 and substituting Equation 13 i ′ ′ = − g i ′ + i s ′ + i ′ s = − g i ′ − i 2 s + s i ′ = − g i ′ + i ′ g + i ′ i − i 2 g + i ′ i = − g i 2 + i ′ ′ 2 i − i i ′

The system admits elementary solutions i = i 0 e − g t , r = r 0 + i 0 1 − e − g t , which correspond to s 0 = 0 . This can be verified by direct substitution into Equation 22. Remarkably, this solution already has the form of a Prony series.

An equivalent second-order system for the s-variable was deduced in [7].

Proposition 2

The SIR system is reducible to the non-linear differential equation for the s-variable s ′ = s ′ 2 s + s − g s ′ or the system S t ′ = S t 2 S + S − g S t (25) S ′ = S t (26)

Proof. d d t i = s − g i = 1 − g s s i = − 1 − g s s ′

On the other hand, d d t i = − d d t s ′ s = s ′ 2 s 2 − s ′ s

Substitution gives s ′ 2 s 2 − s ′ ′ s = − 1 − g s s ′ from where the result follows.

Finally, the r-variable can be represented by the following system.

Proposition 3

The SIR system is reducible to the second order non-linear differential equation for the r-variable r ′ = g r ′ R e e − r g − 1 (27) or the dynamical system R t ′ = g R e R t e − R g − g R t (28) R ′ = R t (29)

Proof. differentiating Equation 15 yields r ′ = g i ′ = g ( s i − g i ) = g s i − g r ′ = − g s ′ − g r ′ . On the other hand, in the r-s plane it holds that d r d s = − g s which possesses an elementary solution as an equation of state: r = − g ⁡ log ⁡ s + c (30)

We use the fixed-point condition s 0 = s ( 0 ) = g so that s = g e − r / g + r 0 / g for the constant r 0 . Further, observe that r 0 = − g ⁡ log g N = g ⁡ log R e where N is the population size from where the result follows.

Obviously, the presented list of dynamical systems Equations 22–29 is non-exhaustive since for any composition of the model variables one could derive a corresponding dynamical system.

The dynamical system described in Proposition 3 is illustrated in Figure 1. The figure exhibits the vector field given by Equations 28, 29 and features one particular trajectory, integrated numerically with the Runge-Kutta fourth order method. The above formulations indicate three potential lines of attack for obtaining a globally convergent series solution. If the system for the r-variable is employed one needs to compute the exponential of a series followed by multiplication with the original series. This seems the least computationally complex approach and it will be pursued in the present paper.

First, we will look for a solution, approximating the r-variable on the positive real line. The asymptotic analysis indicates that an approximation may be achievable using a general Dirichlet series of the form: ρ − ρ ∞ = ∑ k = 1 ∞ c k e − k m τ k ! , τ ≥ 0 (40) where m is a real-valued exponent to be determined later. The truncation of the series at a finite k is then a Prony series. Since we will characterize a dynamical system the Prony series will have exact instead of empirical character. Furthermore, since we are interested only in a Prony approximation we will treat the Equation 40 as formal series without regard to convergence issues.

The formal exponent of the Dirichlet series is another general Dirichlet series e ρ ∞ − ρ = 1 + ∑ k = 1 ∞ − 1 k B k e − k m τ k ! where the coefficients B k are given by the complete, exponential Bell polynomials B k ≡ B k ( c 1 , … , c k ) .

The complete Bell polynomials in turn can be readily calculated from the determinant B k c 1 , … , c k = k − 1 1 c 1 k − 1 2 c 2 k − 1 3 c 3 … c k − 1 k − 2 1 c 1 k − 2 2 c 2 … c n − 1 0 − 1 k − 3 1 c 1 … c k − 2 0 0 − 1 … c k − 3 … … … … … 0 0 0 − 1 c 1

The first few polynomials can be listed as B 0 = 0 , B 1 = c 1 , B 2 = c 2 + c 1 2 , B 3 = c 3 + 3 c 1 c 2 + c 1 3 , B 4 = c 4 + 4 c 1 c 3 + 3 c 2 2 + 6 c 1 2 c 2 + c 1 4 , B 5 = c 5 + 5 c 1 c 4 + 10 c 2 c 3 + 10 c 1 2 c 3 + 15 c 1 c 2 2 + 10 c 1 3 c 2 + c 1 5 , …

We substitute the above expressions in Equation 27, which obtains the form ρ ̈ = ρ ̇ R e e − ρ − 1 Observe that reversing time direction will result in a change of sign of the right-hand side of the equation. Therefore, it will be sufficient to determine the coefficients only for the positive time direction. The resulting expression is m 2 ∑ k = 1 ∞ k 2 c k e − k m τ k ! = − m ∑ k = 1 ∞ k c k e − k m τ k ! R e e − ρ ∞ − 1 + R e e − ρ ∞ ∑ k = 1 ∞ − 1 k B k e − k m τ k !

Since we work with formal infinite series one can apply the Cauchy product formula. Therefore, we obtain the equation R e m ∑ k = 1 ∞ k c k e − ρ ∞ − k m τ k ! + R e m e − ρ ∞ ∑ k = 2 ∞ ∑ i = 1 k − 1 − 1 i B i k − i c k − i e − k m τ i ! k − i ! + m 2 ∑ k = 1 ∞ k 2 c k e − k m τ k ! − m ∑ k = 1 ∞ k c k e − k m τ k ! = 0

Equating the exponents results in the following recursion dependence between the coefficients: k 2 c k m e ρ ∞ − k c k e ρ ∞ + R e ∑ i = 1 k k ! − 1 i B i c k − i i ! k − i − 1 ! + R e k c k = 0 (41)

The first terms of the recursion are 0 , c 1 e ρ ∞ m − e ρ ∞ + R e , 2 2 c 2 e ρ ∞ m − c 2 e ρ ∞ + c 2 R e − c 1 2 R e , 3 3 c 3 e ρ ∞ m − c 3 e ρ ∞ + c 3 R e − c 1 c 2 R e + c 1 3 R e , …

Therefore, for the base exponent m = 1 − R e e − ρ ∞ = 1 + W + − R e e − R e must hold as a constraint for a non trivial series solution to exist. By using the value of m found in this way c k can be obtained by recursion from Equation 41: c k = − q k − 1 ∑ i = 1 k − 1 i k − 1 i B i c k − i , k ≥ 2 where we have recognized the appearance of the Newton’s binomial coefficients and q = R e e ρ ∞ − R e = − W + − R e e − R e 1 + W + − R e e − R e

The first few coefficients of the series can be readily computed as c 1 = λ , c 2 = − q λ 2 , c 3 = q q + 1 λ 3 2 , c 4 = q q 2 − 2 q − 1 λ 4 3 , c 5 = − q q − 1 2 q 2 + 5 q + 1 λ 5 4 , c 6 = − q 27 q 4 − 37 q 3 − 49 q 2 + 49 q + 6 λ 6 30 , c 7 = q 98 q 5 + 311 q 4 − 295 q 3 − 227 q 2 + 197 q + 12 λ 7 72 , …

Substituting the value of q produces an equivalent form of the sequence as c 2 = c 1 2 R e e ρ ∞ − R e , c 3 = − c 1 3 R e e ρ ∞ − 2 R e 2 e ρ ∞ − R e 2 , c 4 = c 1 4 R e e 2 ρ ∞ − 4 R e e ρ ∞ + 2 R e 2 3 e ρ ∞ − R e 3 , c 5 = − c 1 5 R e e ρ ∞ e 2 ρ ∞ − 7 R e e ρ ∞ + 8 R e 2 4 e ρ ∞ − R e 4 , … The above procedure is essentially algorithmic and easily amenable to programming in a computer algebra system. The procedure leaves the coefficient c 1 = λ undetermined. Therefore, the result is a combined series in τ and λ : ρ τ | λ = ρ ∞ + ∑ k = 1 ∞ c k λ k e − k m τ k ! , τ ≥ 0 where now c k are appropriately re-scaled from c 1 = 1 . This amounts to a time-sift of the origin with log ⁡ λ / m since ρ τ | λ = ρ ∞ + ∑ k = 1 ∞ c k e − k m τ − log ⁡ λ k ! , τ ≥ 0 , λ ≠ 0

Furthermore, observe that λ = 0 corresponds to the positive infinity, since then ρ ( τ | 0 ) = ρ ∞ .

In this section we revert the time direction and look for series of the form ρ = ∑ k = 1 ∞ c k e k m τ k ! , τ < 0 which asymptotically approaches ρ − ∞ = 0 . The first terms of the recursion are 0 , c 1 m − R e + 1 , 2 2 c 2 m − c 2 R e + B 1 c 1 R e + c 2 , 3 3 c 3 m − c 3 R e + 2 B 1 c 2 R e − c 1 B 2 R e + c 3 , 4 4 c 4 m − c 4 R e + 3 B 1 c 3 R e + c 1 B 3 R e − 3 B 2 c 2 R e + c 4 , …

Therefore, m = R e − 1 . The same kind of analysis leads to the recursion k 2 c k m − R e k ! ∑ i = 1 k − 1 − 1 i B i k − i c k − i i ! k − i ! − R e k c k + k c k = 0

Therefore, c k = R e R e − 1 k − 1 ∑ i = 1 k − 1 i k − 1 i B i c k − i , k ≥ 2

This results in coefficients c 2 = c 1 2 R e R e − 1 , c 3 = − c 1 3 R e 2 R e − 1 2 R e − 1 2 , c 4 = c 1 4 R e 2 R e 2 − 4 R e + 1 3 R e − 1 3 , c 5 = − c 1 5 R e 8 R e 2 − 7 R e + 1 4 R e − 1 4 , c 6 = c 1 6 R e 4 R e 4 + 36 R e 3 − 134 R e 2 + 73 R e − 6 30 R e − 1 5 , …

Therefore, as before we can set c 1 = λ and obtain the series ρ τ | λ = ∑ k = 1 ∞ c k λ k e k m τ k ! , τ < 0

In summary, as before we can determine the solution up to translation in time.

In this section, we use the peak parametrization where ρ 0 = ρ m = log R e and form the infinite series ρ − ρ 0 ≔ ∑ k = 1 ∞ c k 1 − e − k m τ k ! = A − ∑ k = 1 ∞ c k e − k m τ k ! , t ≥ 0 where we also denote the infinite sum A ≔ ∑ k = 1 ∞ c k / k ! . Then, at positive infinity, the value of A can be determined as ρ ∞ − ρ 0 = R e + W + − R e e − R e − log R e = A ≥ 0

The Lambert function identity W + − 1 / e = − 1 shows that R e = 1 implies A = 0 .

The same development procedure produces the recursion c k = − q k − 1 ∑ i = 1 k − 1 i k − 1 i B i c k − i , k ≥ 2 (42) where now q = 1 e A − 1

The coefficients are still given by Equation 42 but expressed in terms of A can be read off as c 1 = λ , c 2 = λ 2 e A − 1 , c 3 = − λ 3 e A − 2 2 e A − 1 2 , c 4 = λ 4 e 2 A − 4 e A + 2 3 e A − 1 3 , c 5 = − λ 5 e A e 2 A − 7 e A + 8 4 e A − 1 4 , c 6 = λ 6 6 e 4 A − 73 e 3 A + 134 e 2 A − 36 e A − 4 30 e A − 1 5 , … while the base exponent is m = 1 − e − A .

We start from the equation ρ − ρ 0 ≔ − A + ∑ k = 1 ∞ c k e k m τ k ! , τ < 0 (43) where ρ 0 = log R e . Then asymptotically ρ − ∞ − log R e = − A ⇒ A = log R e

So it follows that ρ = ∑ k = 1 ∞ c k e k m τ k ! ∑ k = 1 ∞ λ k c k k ! = log R e where m = R e − 1 , q = R e / ( R e − 1 ) , and c 1 = λ .

In summary, the presented approach is manifestly time-asymmetric due to the properties of the equivalent dynamical system for the r-variable.

The other observable of the model is the ι variable, which can be obtained from Equation 32: ι τ = R e − σ τ − ρ τ = R e 1 − exp − ρ τ − ρ τ

Therefore, by substitution ι τ = R e 1 − exp − ρ 0 + A + ∑ k = 1 ∞ c k e − k m τ k ! − ρ 0 + A + ∑ k = 1 ∞ c k e − k m τ k !

For the value at τ = 0 one obtains ι 0 = R e 1 − exp − ρ 0 − ρ 0 = R e − 1 − log R e = ι m and finally ι τ = R e 1 − e − ρ ∞ ⁡ exp ∑ k = 1 ∞ c k e − k m τ k ! − ρ ∞ + ∑ k = 1 ∞ c k e − k m τ k ! , τ ≥ 0

For the case τ < 0 we proceed in a similar way. ι τ = R e 1 − exp − ∑ k = 1 ∞ c k e k m τ k ! + ∑ k = 1 ∞ c k e k m τ k ! , τ < 0

Then ι 0 = R e 1 − e − log R e − log R e = ι m so that the series agree.

The asymptotes at infinity are verified by direct substitution ι ± ∞ = R e 1 + W ± − R e e − R e R e − R e + W ± − R e e − R e = 0

The series for the s-variable can be determined from the state equations as σ τ = R e ⁡ exp − ρ τ = R e e − ρ ∞ ⁡ exp ∑ k = 1 ∞ c k e − k m τ k ! , τ ≥ 0 and σ τ = R e ⁡ exp − ∑ k = 1 ∞ c k e k m τ k ! , τ < 0

Starting from the 0^th order approximation i ( 0 ) ≈ i 0 , it follows that s ( 0 ) ≈ s 0 e − i 0 τ . However, this does not guarantee convergence of the iteration. To establish convergence we observe that s = g is a fixed point of Equation 22 since d i / d s = 0 for this point and, therefore, the Banach theorem can be applied. Therefore, we must take s 0 = g as an initial condition. This corresponds to the peak-value parameterization so that i 0 = i m , i ′ ( 0 ) = 0 and by Equation 27: i τ = i m exp g ∫ 0 τ e − ∫ 0 z i y d y d z − g τ can be formulated as a functional integral equation to be approximated by DJM.

From there the first order approximation for the i-variable becomes the double exponential function i 1 = i m e g i m 1 − e − i m t − g t (44)

The corresponding solution for r will be r 1 = R e g − i 1 + g W ± − R e e i 1 g − R e (45) where we take the non-principal branch W − for t < 0 and the principal one W + for t > 0 . The s-variable can be expressed accordingly s 1 = − g W ± − R e e i 1 g − R e (46)