Analytical solution of the SIR-model for the temporal evolution of epidemics: part B. Semi-time case

As in part A we assume that the infection and recovery rates have the same time dependence, so that their positive ratio k, the inverse reproduction number, is time-independent and constant:

$\begin{equation}K\left(t\right)=k=\frac{{\mu }_{0}}{{a}_{0}}=\frac{1}{{r}_{0}}\in \left(0,\infty \right).\end{equation} \tag{ 11 }$

As indicated k = 1/r₀ equals the inverse of the SIR-basic reproduction number r₀. We emphasize that this special case includes the standard case used by most analysis before that both, the infection and recovery rates, are constants with respect to time.

The assumption (11) 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 due to assumption (11). Then it is convenient to introduce for arbitrary time dependence of the infection rate a(t) the new time variable

$\begin{equation}\tau \left(t\right)={\int }_{0}^{t}\mathrm{d}\xi \enspace a\left(\xi \right),\end{equation} \tag{ 12 }$

so that

$\begin{equation}\frac{1}{a\left(t\right)}\frac{\mathrm{d}}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}\tau }.\end{equation} \tag{ 13 }$

Equations (4), (5) and (7) then read

$\begin{equation}\frac{\mathrm{d}\mathrm{ln}\enspace I\left(\tau \right)}{\mathrm{d}\tau }=S\left(\tau \right)-k,\end{equation} \tag{ 14 }$

$\begin{align}\hfill I\left(\tau \right)& =-\frac{\mathrm{d}\mathrm{ln}\enspace S\left(\tau \right)}{\mathrm{d}\tau },\hfill \\ \hfill \frac{\mathrm{d}R\left(\tau \right)}{\mathrm{d}\tau }& =kI\left(\tau \right)=-\frac{\mathrm{d}}{\mathrm{d}\tau }\left[k\mathrm{ln}\enspace S\left(\tau \right)\right],\hfill \end{align} \tag{ 15 }$

respectively. In equation (15) (dln S/dτ)_τ=0 = −η in order to meet the initial condition (2) for I(0) = η. Equation (15) can be written as

$\begin{align}\hfill \frac{\mathrm{d}}{\mathrm{d}\tau }\left[R\left(\tau \right)+k\enspace \mathrm{ln}\enspace S\left(\tau \right)\right]=0,\end{align} \tag{ 16 }$

or with the integration constant c₀

$\begin{equation}R\left(\tau \right)={c}_{0}-k\enspace \mathrm{ln}\enspace S\left(\tau \right).\end{equation} \tag{ 17 }$

Because of the initial conditions (2) for S(0) = 1 − η and R(0) = 0 the integration constant equals c₀ = k ln(1 − η), so that

$\begin{align}\hfill R\left(\tau \right)& =k\left[\mathrm{ln}\left(1-\eta \right)-\mathrm{ln}\enspace S\left(\tau \right)\right]=k\enspace \mathrm{ln}\enspace \frac{1-\eta }{S\left(\tau \right)}=-k\left[\varepsilon +\mathrm{ln}\enspace S\left(\tau \right)\right]\hfill \end{align} \tag{ 18 }$

or

$\begin{equation}S\left(\tau \right)=\left(1-\eta \right){\mathrm{e}}^{-\frac{R\left(\tau \right)}{k}}.\end{equation} \tag{ 19 }$

From the invariant $\dot {J}\left(t\right)\mathrm{d}t=j\left(\tau \right)\mathrm{d}\tau$ we obtain with equation (12) in the form dτ/dt = a(t) for differential rate of newly infected persons from the desease

$\begin{equation}j\left(\tau \right)=\frac{\mathrm{d}J\left(\tau \right)}{\mathrm{d}\tau }=\frac{\dot {J}\left(t\right)}{a\left(t\right)}=S\left(\tau \right)I\left(\tau \right)=-\frac{\mathrm{d}S\left(\tau \right)}{\mathrm{d}\tau }\end{equation} \tag{ 20 }$

with the initial value j(0) = η(1 − η). The cumulative distribution corresponding to j(τ) is given by

$\begin{align}\hfill J\left(\tau \right)& ={\int }_{-\infty }^{\tau }\mathrm{d}{\tau }^{\prime }\enspace j\left({\tau }^{\prime }\right)=J\left(0\right)-{\int }_{0}^{\tau }\mathrm{d}{\tau }^{\prime }\frac{\mathrm{d}S\left({\tau }^{\prime }\right)}{\mathrm{d}{\tau }^{\prime }}\hfill \\ \hfill & =J\left(0\right)+S\left(0\right)-S\left(\tau \right)=1-S\left(\tau \right),\hfill \end{align} \tag{ 21 }$

where the initial condition J(0) = I(0) = η had been used, in accord with equation (2).

As in part A we introduce the quantity G(τ) by

$\begin{equation}G\left(\tau \right)=-\mathrm{ln}\enspace S\left(\tau \right),\qquad S\left(\tau \right)={\mathrm{e}}^{-G\left(\tau \right)}\end{equation} \tag{ 22 }$

with the initial condition G(0) = −ln(1 − η) = ɛ. According to equations (15) and (18) then

$\begin{equation}I=\frac{\mathrm{d}G}{\mathrm{d}\tau },\qquad R=kG+k\enspace \mathrm{ln}\left(1-\eta \right)=k\left(G-\varepsilon \right),\end{equation} \tag{ 23 }$

so that the sum constraint (1) provides the dynamical equation

$\begin{equation}\frac{\mathrm{d}G}{\mathrm{d}\tau }=1+k\varepsilon -kG-{\mathrm{e}}^{-G}\end{equation} \tag{ 24 }$

or

$\begin{equation}\frac{\mathrm{d}G/\mathrm{d}\tau }{1+k\varepsilon -{\mathrm{e}}^{-G}-kG}=1.\end{equation} \tag{ 25 }$

Integrating equation (25) over τ then provides

$\begin{equation}\tau +{c}_{1}={\int }_{0}^{G}\frac{\mathrm{d}x}{1+k\varepsilon -{\mathrm{e}}^{-x}-kx},\end{equation} \tag{ 26 }$

where the integration constant c₁ has to be determined from the initial condition G(0) = ɛ, so that

$\begin{equation}{c}_{1}={\int }_{0}^{\varepsilon }\frac{\mathrm{d}x}{1+k\varepsilon -{\mathrm{e}}^{-x}-kx}.\end{equation} \tag{ 27 }$

With equation (27) we obtain for equation (26) the formal solution

$\begin{equation}\tau ={\int }_{\varepsilon }^{G}\frac{\mathrm{d}x}{1+k\varepsilon -{\mathrm{e}}^{-x}-kx}\end{equation} \tag{ 28 }$

with $\tau ={\int }_{0}^{t}\mathrm{d}\xi \enspace a\left(\xi \right)$ according to equation (12).

This solution (28) generalizes the known analytical solutions in the literature [2, 3, 29] as it holds for arbitrary time-dependence of the infection rate a(t). The mentioned known solutions can be reproduced with equation (28) by setting τ = a₀ t on its left-hand side resulting from a constant injection rate a₀. As non-pharmaceutical interventions (NPIs) during the pandemic wave generate a time-changing infection rate, the generalized solution (28) is highly valuable to assess quantitatively the effect of the NPIs on the time evolution of the disease. Substituting y = 1 − e^−x and using equation (21) in the form

$\begin{equation}J=1-{\mathrm{e}}^{-G}\end{equation} \tag{ 29 }$

we obtain for the solution (28)

$\begin{equation}\tau ={\int }_{\eta }^{J}\enspace \frac{\mathrm{d}y}{N\left(y\right)},\end{equation} \tag{ 30 }$

$\begin{equation}N\left(y\right)=\left(1-y\right)\left[y+k\varepsilon +k\enspace \mathrm{ln}\left(1-y\right)\right],\end{equation} \tag{ 31 }$

where we have introduced the inverse integrand N(y), that is parameterized by k and η. Taking the derivative with respect to τ highlights the fact that the inverse integrand in (30) is nothing but the differential rate of newly infected persons in terms of J(τ)

$\begin{align}\hfill j\left(\tau \right)={\left(\frac{\mathrm{d}\tau }{\mathrm{d}J}\right)}^{-1}& =N\left(J\right)=\left(1-J\right)\left[J+k\varepsilon +k\enspace \mathrm{ln}\left(1-J\right)\right].\hfill \end{align} \tag{ 32 }$

Once J(τ) and j(τ) are determined from equations (30) and (32), respectively, one obtains for the other interesting quantities according to equations (21), (20) and (18) that

$\begin{align}\hfill S\left(\tau \right)& =1-J\left(\tau \right),\hfill \\ \hfill I\left(\tau \right)& =J\left(\tau \right)+k\varepsilon +k\enspace \mathrm{ln}\left[1-J\left(\tau \right)\right],\hfill \\ \hfill R\left(\tau \right)& =-k\enspace \mathrm{ln}\left[1-J\left(\tau \right)\right]-k\varepsilon .\hfill \end{align} \tag{ 33 }$

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

$\begin{equation}1+k\varepsilon -{\mathrm{e}}^{-{G}_{\infty }}-k{G}_{\infty }=0.\end{equation} \tag{ 34 }$

Consequently,

$\begin{equation}{G}_{\infty }=\frac{1}{k}+\varepsilon +{W}_{0}\left(\alpha \right),\end{equation} \tag{ 35 }$

where W₀ is the principal solution of Lambert's equation and α given in terms of k and ɛ as

$\begin{equation}\alpha =-\frac{{\mathrm{e}}^{-\frac{1}{k}-\varepsilon }}{k}=-\frac{\left(1-\eta \right){\mathrm{e}}^{-1/k}}{k}.\end{equation} \tag{ 36 }$

The argument α of the Lambert's W₀ function in equation (35) is negative for all k and . Accordingly, W₀ ∈ [−1, 0] is negative as well, while G_∞ ⩾ 0. As ɛ ≪ 1 the maximum G_∞(k) is predominantly determined by the inverse basis reproduction number k. In figure 1 we show the maximum G_∞ as a function of k and three different small values of η ≃ ɛ. It can be seen that G_∞ has values significantly greater than unity for k ≪ 1, while it reaches ${G}_{\infty }\approx \sqrt{2\varepsilon }+O\left(\varepsilon \right)$ at k = 1. 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. Note, that the α defined by equation (36) is just 1 − η = S(0) times the α that occurred in part A. This will have qualitative consequences for the asymptotic behaviors, to be discussed later below.

Zoom In Zoom Out Reset image size

For large values of k ≫ 1 the argument (36) is negative and |α(k ≫ 1)| ⩽ (1 − η)/k ≪ 1. With the asymptotic expansion (A-G10) of the Lambert function for small arguments W₀(|z| ≪ 1) ≃ z we obtain in this limit

$\begin{equation}{G}_{\infty }\left(k\gg 1\right)\simeq \varepsilon +\frac{1}{k}-\frac{{\mathrm{e}}^{-\left(\varepsilon +1k\right)}}{k}\simeq \varepsilon \left(1+\frac{1}{k}+\frac{1}{{k}^{2}}\right),\end{equation} \tag{ 37 }$

in excellent agreement with the numerical results for G_∞ shown in figure 1 for values of k ⩾ 2.

The knowledge of G_∞ from equation (35) and the solutions (28) and (30) allow us to derive a number of exact results. Equation (29) readily provides

$\begin{align}\hfill {J}_{\infty }& =1-{\mathrm{e}}^{-{G}_{\infty }}=1-{\mathrm{e}}^{-\frac{1}{k}-\varepsilon -{W}_{0}\left(\alpha \right)}=1+k{W}_{0}\left(\alpha \right),\hfill \end{align} \tag{ 38 }$

where we used the identity e^−aW(z) = (W(z)/z)^a . Likewise, according to equation (30) J_∞ fulfills the transcendental equation

$\begin{equation}{J}_{\infty }+k\varepsilon +k\enspace \mathrm{ln}\left(1-{J}_{\infty }\right)=0,\end{equation} \tag{ 39 }$

which is consistent with equation (38). Then according to equations (32) and (33)

$\begin{equation}{j}_{\infty }=j\left(\tau =\infty \right)=0,\qquad {I}_{\infty }=I\left(\tau =\infty \right)=0\end{equation} \tag{ 40 }$

and

$\begin{align}\hfill {S}_{\infty }& =S\left(\tau =\infty \right)=1-{J}_{\infty }={\mathrm{e}}^{-{G}_{\infty }}=1-k\left({G}_{\infty }-\varepsilon \right),\hfill \\ \hfill {R}_{\infty }& =R\left(\tau =\infty \right)=1-{I}_{\infty }-{S}_{\infty }=k\left({G}_{\infty }-\varepsilon \right),\hfill \end{align} \tag{ 41 }$

where we have used equation (34). In contrast to the all-time case in part A the final values J_∞, R_∞ and S_∞ depend not only on the inverse basic reproduction number k but also on the fraction of initially infected persons ɛ = −ln(1 − η). Inserting G_∞ from equation (35) we obtain with equation (38)

$\begin{align}\hfill {R}_{\infty }& =1+k{W}_{0}\left(\alpha \right)={J}_{\infty },\hfill \\ \hfill {S}_{\infty }& =1-{R}_{\infty }\left(k\right)=-k{W}_{0}\left(\alpha \right)=1-{J}_{\infty }.\hfill \end{align} \tag{ 42 }$

For large values of k ⩾ 5 we note

$\begin{align}\hfill {J}_{\infty }\left(k\gg 1\right)& ={R}_{\infty }\left(k\gg 1\right)\simeq {G}_{\infty }\left(k\gg 1\right)\simeq \varepsilon \left(1+\frac{1}{k}\right),\hfill \\ \hfill {S}_{\infty }\left(k\gg 1\right)& \simeq 1-\varepsilon \left(1+\frac{1}{k}\right),\hfill \end{align} \tag{ 43 }$

where we used equation (42). In figure 1 we plot the fractions R_∞ and S_∞ from (41) as a function of k ∈ [0, 10].

If the final values were known, one can calculate the corresponding k upon inverting the relationship between S_∞ and k as

$\begin{equation}k=\frac{{S}_{\infty }-1}{\mathrm{ln}\left({S}_{\infty }\right)+\varepsilon }.\end{equation} \tag{ 44 }$

Equation (32) can also be used, even without the explicit inversion of equation (30), to derive generally valid expressions for the time of maximum and the maximum level of the differential rate (20). The maximum of the differential rate (20) occurs when the derivative ${\left(\mathrm{d}j/\mathrm{d}J\right)}_{{J}_{0}}=0$ vanishes. With equation (32) we find

$\begin{equation}\frac{\mathrm{d}j}{\mathrm{d}J}=1-2J-k\left(1+\varepsilon \right)-k\enspace \mathrm{ln}\left(1-J\right),\end{equation} \tag{ 45 }$

yielding for J₀ the transcendental equation

$\begin{equation}-\frac{2{J}_{0}}{k}=\mathrm{ln}\left(1-{J}_{0}\right)+\varepsilon +1-\frac{1}{k},\end{equation} \tag{ 46 }$

or

$\begin{equation}{\mathrm{e}}^{-\frac{2{J}_{0}}{k}}=-{\mathrm{e}}^{\varepsilon +1-\frac{1}{k}}\left[{J}_{0}-1\right].\end{equation} \tag{ 47 }$

With the methods detailed in appendix G of part A the analytic solution of this equation can be expressed in terms of the non-principal solution W₋₁ of Lambert's equation

$\begin{equation}{J}_{0}=1+\frac{k}{2}{W}_{-1}\left({\alpha }_{0}\right),\quad {\alpha }_{0}=\frac{2\alpha }{e}\end{equation} \tag{ 48 }$

with α from (36).

Inserting equation (48) in equation (32) and making use of equation (46) yields for the maximum value in reduced time

$\begin{align}\hfill {j}_{\mathrm{max}}=j\left({J}_{0}\right)& =\left(1-{J}_{0}\right)\left(1-{J}_{0}-k\right)=\frac{{k}^{2}}{4}\left\{{\left[1+{W}_{-1}\left({\alpha }_{0}\right)\right]}^{2}-1\right\}.\hfill \end{align} \tag{ 49 }$

All these expressions for J₀, j_max, and the J_∞ etc are exact, as there are no approximations involved. Since W₋₁(α₀) < − 2 for all k, the j_max is positive. But the time, at which j_max is reached, need not be at positive τ, i.e. it may have passed already. As the SIR model of the present part B is not valid for negative times, we have to simply ignore j_max for such cases, as opposed to part A where j_max was uniquely determined by any choice of initial conditions.

However, for a maximum to occur at finite positive times τ₀ > 0, the derivative dj/dτ has to be positive at times 0 ⩽ τ < τ₀. With equations (45) and (46) we readily find

$\begin{align}\hfill \frac{\mathrm{d}j}{\mathrm{d}\tau }& =\frac{\mathrm{d}J}{\mathrm{d}\tau }\frac{\mathrm{d}j}{\mathrm{d}J}=j\frac{\mathrm{d}j}{\mathrm{d}J}=j\left[1-2J-k\left(1+\varepsilon \right)-k\enspace \mathrm{ln}\left(1-J\right)\right]\hfill \\ \hfill & =j\left[2\left({J}_{0}-J\right)-k\enspace \mathrm{ln}\frac{1-J}{1-{J}_{0}}\right].\hfill \end{align} \tag{ 50 }$

Since the requirement of a maximum j_max to exist at positive times is identical to the requirement of a positive dj/dτ at τ = 0, we can insert J(0) = η into the second line of equation (50) to find

$\begin{equation}k{< }{\left.\frac{1-2J}{1+\varepsilon +\mathrm{ln}\left(1-J\right)}\right\vert }_{\tau =0}=1-2\eta ,\end{equation} \tag{ 51 }$

where the relationship (3) between η and ɛ had been used. For inverse reproduction numbers k greater than 1 − 2η, the daily rate is monotonically decreasing at all times from its initial value j(0) = η(1 − η). Consequently, depending on the values of k and η we have two different cases:

• (a)  
  The 'decay' case for large values of k ⩾ 1 − 2η, where the daily rate of newly injected persons monotonically decreases at all positive times from its initial maximum value j(0),
• (b)  
  The 'peak case' if k < 1 − 2η where the daily rate of newly infected persons attains a maximum at a finite positive time.In table 1 we summarize the main relations between the various epidemical quantities as well as their minimum and maximum values.

Table 1. An attempt to help the reader to guide through this document. Peak time is meant to be the time where j reaches its maximum. Note that J₀, equation (48), is very well approximated also by J_*, equation (62). The coefficients a_0,1,2 entering the quoted equations are given by equation (58), coefficients b_0,1,2 are given by equation (63) (version I) and equation (65) (version II), remaining constants written down in equation (72). The exact τ in terms of J is given by equation (30). The peak time does not exist in the future when k < 1 − 2η, according to equation (51), or equivalently, when τ_* < 0. The two versions I and II distinguish themselves only during decay times, and differ qualitatively mainly in their asymptotic behaviors at large times τ ≫ τ_*. The reduced time τ is given in terms of physical time t via equation (12), which allows to model arbitrary time-dependent infection rates a(t).

Quantity X =	S	I	R	J	j	G
In terms of J	1 − J	Eq. (33)	Eq. (33)	id.	Eq. (32)	Eq. (22)
Initial value X(0)	1 − η	η	0	η	η(1 − η)	ɛ
Minimum value Xmin	S∞	η	0	η	0	ɛ
Maximum value Xmax	1 − η	1 − k[1 −ɛ + ln k]	R∞	Eq. (41)	Eq. (49)	G∞
Final value X∞	Eq. (41)	0	Eq. (41)	Eq. (38)	0	Eq. (35)
Peak time τmax ≃ τ*	—		—	—	Eq. (73)	—
Value at peak time	1 − J0	Eq. (33)	Eq. (33)	J0	jmax	Eq. (22)
In terms of τ (τ ⩽ τ*)	Via J	Via J	Via J	Eq. (70)	Eq. (74)	Via J
In terms of τ (τ ⩾ τ*)	Via J	Via J	Via J	Eq. (71)	Eq. (74)	Via J
Asymptotic τ → ∞ behavior					${\mathrm{e}}^{\left(1-k-{J}_{\infty }\right)\tau }$

An ideal approximant of the integrand must have several features:

• (a)  
  The integral τ(J) must exist in analytic form,
• (b)  
  The τ(J) function must exhibit an analytic inverse J(τ),
• (c)  
  The exact and approximated reciprocal integrands must coincide at y = η to make sure the initial j(τ) is correctly captured,
• (d)  
  The exact and approximated reciprocal integrands must coincide at y = J_∞ to make sure J approaches the exact J_∞ at large times,
• (e)  
  The exact and approximated reciprocal integrands should have the same slope in the vicinity of the integral bounds, as these coefficients determine the asymptotic behavior of all SIR functions.If the integrand is approximated by two functions over disjunct regimes (a) and (b) separated by some y value (denoted by J_*) there are further requirements:
• (f)  
  The two approximants must possess identical values at y = J_* to prevent a discontinuity in J(τ), and
• (g)  
  They should have also the same first derivative at y = J_* to prevent a discontinuity in j(τ) when calculated via J'(τ).Requirement (g) can be circumvented upon using the exact relation (32) or alternatively, and more conveniently, by just using j(τ) as given by the approximant, i.e., j(τ) = j_a (J(τ)) within regime (a), and j(τ) = j_b (J(τ)) within regime (b). We shall present below two approximants that both offer all but one of these ideal features. While the first approximant (referred to a version I) has (a)–(f), the second approximant (referred to as version II) has (a)–(d) and (f) and (g). The situation is depicted for both versions in figure 2.

Zoom In Zoom Out Reset image size

To account for (a) and (b), with a side-view to the remaining requests, we note that the following indefinite integral, carrying a 2nd order polynomial in the denominator of the integrand (with the abbreviation ${c}_{3}=\sqrt{{c}_{1}^{2}-4{c}_{0}{c}_{2}}$),

$\begin{align}\hfill {T}_{c}\left(J\right)& \equiv {\int }^{J}\frac{\mathrm{d}y}{{\sum }_{j=0}^{2}{c}_{j}{\left(y-{y}_{c}\right)}^{j}}=-\frac{2}{{c}_{3}}\enspace {\mathrm{tanh}}^{-1}\left[\frac{{c}_{1}+2{c}_{2}\left(J-{y}_{c}\right)}{{c}_{3}}\right],\hfill \end{align} \tag{ 54 }$

can be inverted analytically to yield

$\begin{equation}{J}_{c}\left(T\right)={y}_{c}-\frac{{c}_{1}+{c}_{3}\enspace \mathrm{tanh}\left({c}_{3}T/2\right)}{2{c}_{2}}.\end{equation} \tag{ 55 }$

The derivative is thus

$\begin{equation}{j}_{c}\left(T\right)=\frac{\mathrm{d}{J}_{c}\left(T\right)}{\mathrm{d}T}=-\frac{{c}_{3}^{2}}{4{c}_{2}\enspace {\mathrm{cosh}}^{2}\left({c}_{3}T/2\right)}.\end{equation} \tag{ 56 }$

We will make use of these identities in the following section. Here, c is a placeholder for either (a) rise or (b) decay.

To account for the requirements (c) and (e) we use an optimized Taylor-expansion, i.e., Taylor expansion subject to a constraint [33], of the reciprocal integrand in equation (53) about y = y_a = η. Up to 3rd order in y − η one has

$\begin{align}\hfill N\left(y\right)& ={j}_{a}\left(y\right)+\mathcal{O}{\left(y-{y}_{a}\right)}^{3}\hfill \\ \hfill {j}_{a}\left(y\right)& =\sum\limits _{i=0}^{2}{a}_{i}{\left(y-\eta \right)}^{i},\hfill \end{align} \tag{ 57 }$

where the coefficients follow from these requirements

$\begin{equation}{a}_{0}=\left(1-\eta \right)\eta ,\end{equation} \tag{ 58 }$

$\begin{equation}{a}_{1}=1-k-2\eta ,\end{equation} \tag{ 59 }$

$\begin{equation}{a}_{2}=\frac{{j}_{\mathrm{max}}-{a}_{0}-{a}_{1}\left({J}_{0}-\eta \right)}{{\left({J}_{0}-\eta \right)}^{2}}.\end{equation} \tag{ 60 }$

While a₀ and a₁ arise from the Taylor expansion, a₂ ensures that equation (57) evaluated at y = J₀ yields j_max, where both J₀ and j_max are known from equations (48) and (49). This a₂ is only slightly different from the a₂ one obtains via the Taylor expansion. While a₀ is positive, a₂ is negative, and a₁ may have either sign, depending on k. The abbreviation ${a}_{3}=\sqrt{{a}_{1}^{2}-4{a}_{0}{a}_{2}}$ will be used.

To account for (d) and (e) we also expand the reciprocal integrand about y = y_b = J_∞ to write, again up to 3rd order,

$\begin{align}\hfill N\left(y\right)& ={j}_{b}\left(y\right)+\mathcal{O}{\left(y-{y}_{b}\right)}^{3}\hfill \\ \hfill {j}_{b}\left(y\right)& =\sum\limits _{i=1}^{2}{b}_{i}{\left(y-{J}_{\infty }\right)}^{i}\hfill \end{align}, \tag{ 61 }$

where the sum does not involve b₀ as the denominator N(y) vanishes at y = J_∞. For the same reason b₃ = |b₁|. We could mention here the exact values for b₁ and b₂ resulting from Taylor expansion. But this would not help solving our problem. If we would go for using the approximant j_a (y) within some range of y ∈ [η, J_*] values, and j_b (y) within the remaining range y ∈ [J_*, J_∞], the approximant would not share values and derivatives at y = J_*, if the Taylor coefficients are used, and thus violate requirements (f) and (g).

Our solution to this discontinuity problem stems from the observation that j_a (y) with coefficients given by equation (58) approximates the reciprocal integrand very well up to its extremum j_a (J_*) located at

$\begin{equation}{J}_{{\ast}}=\eta -\frac{{a}_{1}}{2{a}_{2}}.\end{equation} \tag{ 62 }$

In accord with our previous statements, a peak time occurs at positive times only when J_* > η implying that a₁ has to be positive (because a₂ is negative), corresponding to k < 1 − 2η in light of (59). We note that this requirement agrees exactly with our earlier condition (51).

To fulfill condition (f) we need to require that j_b (J_*) = j_a (J_*). If we ignore the requirement (g) but take into account (e), then we arrive at our version I,

$\begin{equation}{b}_{1}=1-k-{J}_{\infty },\end{equation} \tag{ 63 }$

$\begin{equation}{b}_{2}=\frac{{a}_{2}\left[2{a}_{1}{b}_{1}-{a}_{3}^{2}+4{a}_{2}{b}_{1}\left({J}_{\infty }-\eta \right)\right]}{{\left[{a}_{1}+2{a}_{2}\left({J}_{\infty }-\eta \right)\right]}^{2}}\end{equation} \tag{ 64 }$

where b₁ is a true Taylor coefficient. The abbreviation ${a}_{3}=\sqrt{{a}_{1}^{2}-4{a}_{0}{a}_{2}}$ has been used. Note that b₁ ⩽ 0 while b₂ may have either sign, depending on k and η. This version thus implies correct asymptotic behavior at large y at the expense of a jump in the derivative of J with respect to τ upon approaching J_* from different sides, J > J_* and J < J_*. This jump does however not occur here, as we do not use a derivative to calculate j from J at J = J_*, but two functions that meet at J_*; the continuous j is already given by the approximants of the inverse integrand, j_a and j_b , which do not possess any jumps, and meet exactly where we split the integral, as this was the basic requirement leading to the coefficients in the present paragraph. Because this version I exhibits an exact Taylor coefficient b₁, that determines the behavior in the vicinity of J_∞, or equivalently, in the long time limit τ → ∞, we can deduce the asymptotic behavior of the exact solution of the semi-time SIR model in section 5.4.

The opposite case of fulfilling (g) while ignoring (e) yields version II

$\begin{equation}{b}_{1}=\frac{{a}_{3}^{2}}{{a}_{1}+2{a}_{2}\left({J}_{\infty }-\eta \right)},\end{equation} \tag{ 65 }$

$\begin{equation}{b}_{2}=\frac{{a}_{2}{a}_{3}^{2}}{{\left[{a}_{1}+2{a}_{2}\left({J}_{\infty }-\eta \right)\right]}^{2}}.\end{equation} \tag{ 66 }$

Also for this version, b₁ ⩽ 0 and b₂ may have either sign, depending on k and η. This version implies a continuously differentiable J, serves to reproduce the exact J up to times well beyond the peak time, but does not exhibit the correct asymptotic behavior at long times. This formal drawback is something that does not matter as long as k is not very small. We moreover know the asymptotic behavior $j\left(\tau \right)\sim {\mathrm{e}}^{{b}_{1}\tau }$ of the exact solution, that helps to rate when version II begins to fail.

Based on the considerations in the last section we are in the position to immediately write down the approximants for both versions. An approximant consists of two parts, J_rise(τ) to be used at τ ⩽ τ_*, and J_decay(τ) to be used at τ > τ_*, where

$\begin{equation}{\tau }_{{\ast}}=\tau \left({J}_{{\ast}}\right)={T}_{a}\left({J}_{{\ast}}\right)-{T}_{a}\left({y}_{a}\right)\end{equation} \tag{ 67 }$

according to equation (54). We will simplify this expression for τ_* further below. To understand how we write down the approximant using the above equations, we recall that our approximation evaluates the integral τ in equation (52) via

$\begin{equation}\tau =\begin{cases}_{\eta }^{J}\mathrm{d}y/{j}_{a}\left(y\right),\quad \hfill & J{\leqslant}{J}_{{\ast}}\hfill \\ {\int }_{\eta }^{{J}_{{\ast}}}\mathrm{d}y/{j}_{a}\left(y\right)+{\int }_{{J}_{{\ast}}}^{J}\mathrm{d}y/{j}_{b}\left(y\right),\quad \hfill & J{\geqslant}{J}_{{\ast}}\hfill \end{cases},\end{equation} \tag{ 68 }$

where the 2nd order polynomials j_a (y) and j_b (y) are characterized by coefficients a_0,1,2 and b_1,2, respectively, and where all these five coefficients as well as the J_* were chosen to fulfill our criteria for an approximant, and explicitly expressed already in terms of k and η. Using equation (54) and τ_* = τ(J_*) we can write equation (68) as

$\begin{equation}\tau =\begin{cases}_{a}\left(J\right)-{T}_{a}\left(\eta \right),\quad \hfill & 0{\leqslant}\tau {\leqslant}{\tau }_{{\ast}}\hfill \\ {\tau }_{{\ast}}+{T}_{b}\left(J\right)-{T}_{b}\left({J}_{{\ast}}\right),\quad \hfill & \tau {\geqslant}{\tau }_{{\ast}}\hfill \end{cases}\end{equation} \tag{ 69 }$

because c → a (and y_a = η) has to be used for the rising J < J_* (or equivalently, τ < τ_*) regime, and c → b (and y_b = J_∞) afterwards. For the rising regime τ < τ_* we can thus write τ + T_a (η) = T_a (J), where T_a (η) is a constant. According to equation (55) the inversion is

$\begin{equation}{J}_{\text{rise}}\left(\tau \right)=\eta -\frac{{a}_{1}+{a}_{3}\enspace \mathrm{tanh}\left\{\frac{{a}_{3}}{2}\left[\tau +{T}_{a}\left({y}_{a}\right)\right]\right\}}{2{a}_{2}},\end{equation} \tag{ 70 }$

where y_a = η has been used, and because the τ + T_a (η) has overtaken the role of T in equation (55).

Similarly, within the decay regime τ ⩾ τ_*, equation (69) can be written as τ − τ_* + T_b (J_*) = T_b (J), where τ_* and T_b (J_*) are constants. According to equation (55) the inversion is

$\begin{equation}{J}_{\text{decay}}\left(\tau \right)={J}_{\infty }-\frac{{b}_{1}}{2{b}_{2}}\left[1+\mathrm{tanh}\left\{\frac{{b}_{1}}{2}\left[\tau -{\tau }_{{\ast}}+{T}_{b}\left({J}_{{\ast}}\right)\right]\right\}\right],\end{equation} \tag{ 71 }$

where y_b = J_∞ has been used, and because the τ − τ_* + T_b (J_*) has overtaken the role of T in equation (55).

It seems useful to explicitly write down the constants appearing in the above final expressions

$\begin{align}\hfill {T}_{a}\left({y}_{a}\right)& =-\frac{2}{{a}_{3}}\enspace {\mathrm{tanh}}^{-1}\left(\frac{{a}_{1}}{{a}_{3}}\right),\hfill \\ \hfill {T}_{a}\left({J}_{{\ast}}\right)& =-\frac{2}{{a}_{3}}\enspace {\mathrm{tanh}}^{-1}\left[\frac{{a}_{1}+2{a}_{2}\left({J}_{{\ast}}-\eta \right)}{{a}_{3}}\right]=0\hfill \\ \hfill {T}_{b}\left({J}_{{\ast}}\right)& =-\frac{2}{{b}_{1}}\enspace {\mathrm{tanh}}^{-1}\left[\frac{{b}_{1}+2{b}_{2}\left({J}_{{\ast}}-{J}_{\infty }\right)}{{b}_{1}}\right],\hfill \end{align} \tag{ 72 }$

where we made use of equation (62) to simplify T_a (J_*), and where J_* was already expressed in terms of the coefficients in equation (62). As an aside we emphasize that by construction, more precisely, requirement (f), the two expressions J_rise(τ_*) = J_decay(τ_*) agree. Inserting the coefficients from equation (72) we obtain for the maximum relative time (67) the simplified, still identical expression

$\begin{equation}{\tau }_{{\ast}}=\frac{2}{{a}_{3}}\enspace {\mathrm{tanh}}^{-1}\left(\frac{{a}_{1}}{{a}_{3}}\right)\end{equation} \tag{ 73 }$

so that T_a (y_a ) in equations (72) and (70) turn out to be just identical to −τ_*. These are final expressions for J(τ) in terms of the known coefficients, valid for both versions. If τ_* lies in the past, or equivalently, if k < 1 − 2η, then J(τ) ≃ J_decay(τ) is the only contribution to our approximant.

With this approximant for J(τ) at hand, the differential number of newly infected persons is simply given by j_rise(τ) = j_a [J_rise(τ)] with the 2nd order polynomial j_a from equation (57), and j_decay(τ) = j_b [J_decay(τ)] with the 2nd order polynomial j_b from equation (61). More explicitly,

$\begin{equation}{j}_{\text{rise}}\left(\tau {\leqslant}{\tau }_{{\ast}}\right)=\sum\limits _{i=0}^{2}{a}_{i}{\left[{J}_{\text{rise}}\left(\tau \right)-\eta \right]}^{i},\end{equation} \tag{ 74 }$

$\begin{equation}{j}_{\text{decay}}\left(\tau {\geqslant}{\tau }_{{\ast}}\right)=\sum\limits _{i=1}^{2}{b}_{i}{\left[{J}_{\text{decay}}\left(\tau \right)-{J}_{\infty }\right]}^{i}.\quad \end{equation} \tag{ 75 }$

Because version II has the property (g), the differential number of cases is alternatively given by equation (56), i.e.,

$\begin{equation}{j}_{\text{rise}}\left(\tau \right)=-\frac{{a}_{3}^{2}}{4{a}_{2}\enspace {\mathrm{cosh}}^{2}\left[{a}_{3}\left(\tau -{\tau }_{{\ast}}\right)/2\right]},\end{equation} \tag{ 76 }$

$\begin{equation}{j}_{\text{decay}}\left(\tau \right)=-\frac{{b}_{1}^{2}}{4{b}_{2}\enspace {\mathrm{cosh}}^{2}\left\{{b}_{1}\left[\tau -{\tau }_{{\ast}}+{T}_{b}\left({J}_{{\ast}}\right)\right]/2\right\}}.\end{equation} \tag{ 77 }$

For version I, the daily new number of cases is not properly described by dJ/dτ, but instead is obtained from J via equation (74). For version II, equations (76) and (77) are just identical to equation (32) and also identical to equation (74) when J(τ) from equations (70) and (71) is used.

We recall that the remaining SIR quantities are directly obtained as well from J(τ) via S(τ) = 1 − J(τ), I(τ) = J(τ) + kɛ + k ln[1 − J(τ)], and R(τ) = 1 − S(τ) − I(τ).

While the value j_max of the maximum of j(τ) at peak time τ_max is captured exactly by construction of our approximants, the approximate peak time τ_* given by equation (73) is compared with the numerically exact peak time τ_max of the differential number of cases j(τ) in figure 3(a) for some selected values of η, to appreciate the dependency on k. The relative error for the selected parameter η values is well below 2% for all k. Such relative error of our approximant τ_*, now as function of k and η, i.e., capturing the whole parameter space, is shown in figure 3(c). For these calculations and comparisons, the exact peak time τ_max is obtained by solving equation (30), or alternatively, the SIR time evolution equations (14) and (15) numerically. As is visible from figure 3(c), the relative error of τ_* does not exceed 2.5% in the worst case, which is located close to k = 0.6 and η = 10^−1.5 ≈ 0.03. For both small and large k the relative error is extremely small, well below 0.1%. Our simple analytical approximant τ_* for the peak time of daily newly infected persons can thus be considered precise for all practical purposes. As is obvious from figure 3(a), the peak time decreases with increasing η = I(0), for a given k.

Zoom In Zoom Out Reset image size

The same holds true for the approximant J_* of the cumulative number of infected persons at peak time, for which we know the exact analytic expression J₀ as well. The two are compared with each other in figure 3(b) for some selected values of η, to appreciate the dependency on k. The relative error of our approximant J_* over the whole parameter space is shown in figure 3(d). Note that J_* ⩾ η for the semi-time SIR, while J₀ < η, or equivalently, k > 1 − 2η, marks the regime where the peak time does not lie in the future. Under such circumstances, there is only a regime of decay (J_* = η). Note also, that J₀ can get negative here, as opposed to part A, which corresponds to a choice of initial conditions that are incompatible with the all-time SIR model. More specifically, J₀ < 0 for k < 1 − 2η, which is just fine and irrelevant, while J_* = η within this regime, which is characterized by a purely decaying J(τ). As shown by figure 3(b), the cumulative number of infected persons at peak time decreases with increasing η = I(0), for a given k.

By our construction, the j_rise(τ) achieves its maximum at time τ_*, where j_rise(τ_*) = j_max is identical to the exact analytic result (49). The peak time τ_max of the unapproximated j(τ) (if a peak time exists) is not exactly located at τ_*, because J_* is not exactly identical to J₀ (figure 3(b)), as we discussed already. When there is no maximum in j at positive times, J_* = η.

In figures 4(a) and (b) we plot the daily rate j(τ) of newly infected persons for the more relevant cases where a peak time is ahead in the future. These figures compare the time evolution of our approximants (version I and II) with the exact numerical solution of the SIR model. Decay cases are shown in figures 4(c) and (d), again for selected values of the parameters k and η. To make sure that we did not just pick suitable parameters here, the mean relative error of our approximant for j(τ) is evaluated as function of k and η in figures 4(e) and (f) for both versions, and over the whole range of parameter space. The relative errors of the approximant is below 0.3% everywhere for version II, and below 1% for all k > 0.02 for version I. As already known from figure 3(a) the peak time increases with increasing η and increasing k before dropping sharply when k approaches unity. There is no future peak time anymore for k > 1 − 2η, reflected by a negative τ_* in that case.

Zoom In Zoom Out Reset image size

Similarly, in figure 5 we compare the approximations (70) and (71) for the cumulative number with the exact time dependence J(τ) obtained by the numerical integration of equation (30), for both versions I and II. While figures 5(a) and (b) focuses on the case where the peak times lies in the future, figures 5(c) and (d) is concerned with the decay cases. It is visible how the cumulative number at peak time, J₀, approximated by J_*, decreases with increasing k, in accord with figure 3(b), and is only weakly affected by the initial condition for η ≪ 1. The accuracy of the approximant is below 0.3% for any choice of parameters, if version II is used, while it exceeds 1% for version I but only for small k < 0.1.

Zoom In Zoom Out Reset image size

Relative errors for all other SIR quantities are similarly small, and need not be discussed here, as all the S, I, R, G functions can be derived from J(τ), as table 1 mentioned already.

Overall, the difference between versions I and II is not easily appreciated from these figures, and therefore of minor relevance in practice. Both versions are in any case formally very similar as they only differ in the values of their polynomial coefficients b₁ and b₂, and thus in their behavior beyond the peak time. At early times the doubling time τ₂ defined by j(τ + τ₂) = 2j(τ) is given by τ₂ = ln(2)/a₃, according to equation (76). The difference between versions I and II is most pronounced if one focuses solely on j(τ) at very late times.

The differential rate of newly infected persons j(τ) exhibits an exponential decay

$\begin{equation}\underset{\tau \to \infty }{\mathrm{lim}}j\left(\tau \right)\propto {\text{e}}^{-\nu \tau }\end{equation} \tag{ 78 }$

with some semipositive and constant, dimensionless decay rate ν at late times τ ≫ τ_max under all conditions, that is obtained from equation (63). The decay rate is identical to ln(2)/τ_1/2, if τ_1/2 denotes the late half time, i.e., the time required for j to have dropped by a factor 2. Figure 6 shows the exact rate

$\begin{equation}\nu =k+{J}_{\infty }-1=k\left[1+{W}_{0}\left(\alpha \right)\right]\end{equation} \tag{ 79 }$

of the semi-time SIR as function of k, for various η, in perfect agreement with the decay rate featured by our approximant version I (ν = −b₁ in that case). In equation (79) we have used equation (38) to replace J_∞ by Lambert's function, with α defined by equation (36). This way the decay rate is formally identical to part A, recalling, that α for the present semi-time SIR model is S(0) times the α for the all-time SIR. The rate ν is positive, except for k = 0, and for k = 1 if η = 0. The decay rate generally increases with increasing initial fraction of infected persons, η. This means that the daily rate at late times decreases the faster the larger the initial fraction of infected persons is. Shown for comparison (dashed black) is the corresponding decay rate for the all-time SIR model from part A. We thus find that the exponential decrease is under all conditions faster for the all-time SIR, compared with the semi-time SIR. While for small η the exponent ν is almost symmetric about k = 1, it gets more and more asymmetric upon increasing η. We recall that the all-time SIR model is restricted to the regime k ∈ [0, k_max] with k_max = [S(0) − 1]/ln(S(0)) which is ${k}_{\mathrm{max}}=\eta /\varepsilon \in \left(\right.\enspace 0,1\left. \right]$ using our current notation.

Zoom In Zoom Out Reset image size