Explicit formulae for the peak time of an epidemic from the SIR model. Which approximant to use?

An analytic evaluation of the peak time of a disease allows for the installment of effective epidemic precautions. Recently, an explicit analytic, approximate expression (MT) for the peak time of the fraction of infected persons during an outbreak within the susceptible–infectious–recovered/removed (SIR) model had been presented and discussed (Turkyilmazoglu, 2021). There are three existing approximate solutions (SK-I, SK-II, and CG) of the semi-time SIR model in its reduced formulation that allow one to come up with different explicit expressions for the peak time of the infected compartment (Schlickeiser and Kröger, 2021; Carvalho and Gonçalves, 2021). Here we compare the four expressions for any choice of SIR model parameters and find that SK-I, SK-II and CG are more accurate than MT as long as the amount of population to which the SIR model is applied exceeds hundred by far (countries, ss, cities). For small populations with less than hundreds of individuals (families, small towns), however, the approximant MT outperforms the other approximants. To be able to compare the various approaches, we clarify the equivalence between the four-parametric dimensional SIR equations and their two-dimensional dimensionless analogue. Using Covid-19 data from various countries and sources we identify the relevant regime within the parameter space of the SIR model.


Introduction
The temporal evolution of COVID-19 (or SARS-CoV-2) pandemic waves has been successfully described, discussed, and forecasted by the mathematical susceptible-infectious-recovered/ removed SIR model [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] while the model itself had been developed nearly a century ago [18][19][20][21][22]. There are two noticeable quantities that indicate the occurrence of new pandemic waves: the fraction i(t) of infected persons at time t, and the fraction of newly infected population per dayJ(t) = β(t)s(t)i(t), where β(t) denotes the infection rate, and s(t) the fraction of susceptible population. Both indicators i(t) andJ(t) in a pandemic wave first increase with time, undergo a maximum and drop at late times. The two peak times are slightly different. While the peak time inJ(t) is the one usually reported in the media on the basis of reported number of newly infected persons, the peak time in i(t) is the one that determines the peak time of required clinical resources. While an analytic approximant for the usually measurable peak time inJ exists for both the all-time [18,23] approximants supersedes the one presented to them recently by Turkyilmazoglu [25], or not.
There are various papers [23][24][25][26][27][28][29][30][31][32][33] that aimed at deriving approximate analytic expressions for various quantities that appear in the all-time or semi-time SIR models [23][24][25][26][27][28][29][30][31][32]34,35]. And there are several variants [36][37][38][39][40][41][42][43] of the SIR model, including stochastic variants [44][45][46][47][48][49][50][51][52][53][54][55][56][57] or variants that account for vaccination [58][59][60][61][62] or pulse vaccination [63][64][65][66]. For the semi-time SIR Heng and Althaus [26] provided an analytic approximant for the population fraction of susceptible s(t) and infected i(t) persons but only for times prior peak time, and they could not derive an approximant for any of the peak times. Harko et al. [67] were able to express time of the semi-time SIR model in terms of an integral, which can only be evaluated numerically, but did not derive an integral or analytic expressions for peak times. Bidari et al. [27] focused on deriving implicit equations for final sizes of compartments and approximate solutions to such equations, but did not obtain an expression for a peak time. An analytic approximate inverse solution of the semi-time SIR model with special initial condition of initially vanishing population fraction r(0) = 0 of recovered/removed persons was presented by Carvalho and Gonçalves [28]. An explicit series solution of SIR and SIS epidemic models was obtained by making use of the homotopy analysis method by Khan et al. [30]. This approach allows one to study convergence properties of the solution, but in practice, it is inefficient to study the solution of the SIR model using an infinite series. Moreover, this approach does not provide an analytic expression for a peak time. Within the same spirit Barlow and Weinstein [33] derived an accurate closed-form infinite series solution of the SIR model that involves a Vandermonde matrix, whose inversion is explicitly known; however, they do not provide a method for estimating peak time.
In [25] the time t-dependent SIR model for population fractions s (susceptible), i (infected), and r (removed/recovered) had been written as follows subject to the general arbitrary (semi-time) initial conditions and where β and γ are the assumed stationary infection and recovery rates of population fractions. The initial condition for r follows from r(0) = r 0 = 1 − s 0 − i 0 , because the classical SIR model has only three compartments, and any of the N members of the population belongs to one of the three compartments at any time. Allowing for arbitrary initial conditions and restricting the model to future times t > 0 is known as semi-time SIR model [19][20][21]24]. For the all-time SIR model, s 0 and i 0 are interrelated [18,23]. Here we consider the semi-time SIR model for which an analytic expression for the peak time t peak of the infected compartment, defined by had been derived in [25]. It will be reproduced in Eq. (14). In its classical form Eq. (1) with (2) the dimensional SIR model has four parameters, the rates β and γ , and the initial fractions s 0 and i 0 .
On the other hand, analytic approximants for the solution of the reduced dimensionless semi-time SIR model are available from [24] and [28]. While the reduced model has only three parameters, it is still equivalent with the original SIR model (the proof of equivalence will be given in the next section). Moreover, one of the three parameters is used to make the time dimensionless, so that one is left with essentially only two parameters k and η. These parameters appear in the reduced semi-time SIR model [24] as follows with initial conditions Here in general a dimensionless reduced time τ defined by τ ∝ ∫ t 0 β(ξ )dξ for an arbitrary (including periodic) time-dependent infection rate β(t), but because as in [25] we consider here a stationary infection rate, τ is simply proportional to time, τ = at with a constant rate a that is specified in Eq. (7). Moreover, S = s/(s 0 + i 0 ), I = i/(s 0 + i 0 ), and R = 1 − S − I are fractions of the initially unrecovered population. This implies R(0) = 0 as in [24]. We recall, that the quantities s, i, and r are fractions with respect to the total initial population, including the recovered compartment, and therefore s ≥ S, i ≥ I, and r ≥ R.
Yet another version of a reduced semi-time SIR model was investigated in [28]. They introducedτ = γ t to write subject to unchanged initial conditions (5). It is important to note here that in the original work Carvalho and Gonçalves [28] formulated their equations in terms of the fractions s, i, and r, did not allow for non-vanishing r(0), and mentioned R 0 = β/γ , which is a correct assignment but only for this special case of r(0) = 0.
We are here extending their work to allow for non-vanishing r(0) so that also R 0 has to be revised, as shown next. This extension will allow one to compare all available approximants for arbitrary initial conditions s 0 and i 0 .

Analytical approximants
Starting from the dimensional SIR model with four parameters, the three parameters of the reduced formulation are given by In turn, for given k, η, and a, the parameters of the dimensional model are given upon direct inversion of Eqs. (7) just highlighting the fact, that s 0 must drop out. Moreover, any characteristic real time t of the dimensional SIR model is related to the corresponding reduced time τ = at or alsoτ = γ t as where the dimensionless times τ peak (k, η) andτ peak (R 0 , η) depend on two parameters only: the inverse basic reproduction number k = 1/R 0 and the initial fraction of infected population among the non-recovered population, η = I(0).
If we wish to test an approximate solution of the SIR model for arbitrary choices of parameters, or if we want to compare the quality of different approximants as we are going to do here, it is hence sufficient to perform the test in the 2-dimensional parameter space built by k and η. While η ∈ (0, 1) by construction, the k is semipositive in general. However, a peak time in I(t) occurs at positive times t > 0 only if the following inequality holds This inequality follows from [24], as the condition (3) for the peak converts with the help of (4) into S(τ peak ) = k (11) within the reduced model, because S(0) = 1 − η, and because S monotonically decreases with increasing τ according to Eq. (4).
We have to still prove the equivalence between dimensional and dimensionless forms of the SIR model. To this end it is sufficient to prove Eq. (7), or the equivalent Eq.
confirming the equation of change for I in Eq. (4), which is obvious, if we divide both sides of the Eqs. (12) and (13) by (s 0 + i 0 )a. The initial conditions are also equivalent, as s 0 = (s 0

MT approximant
The MT approximant by Turkyilmazoglu [25] for the peak time of the infected compartment had been formulated using the dimensional SIR formulation (1) with 4 parameters. Using our replacement rule (8), it receives the form of Eq. (9) with the dimensionless peak time (see Appendix A for details) where k-and η-dependent coefficients are given by Note that the term under the square root E is positive for all η ∈ (0, 1) and k ∈ (0, 1). For the special case of k = 2/3, where the denominator of τ MT peak vanishes, the expression can still be evaluated using l'Hopital's rule, i.e., τ MT

SK approximants SK-I and SK-II
The Schlickeiser & Kröger (SK) approximants for the peak time of the infected compartment, starting from the reduced SIR model Eqs. (4) had not been explicitly written down in their work [24], where they developed an approximate solution to the whole time-dependency of all SIR quantities, but they can be read off from the provided approximate analytic solution S(τ ) of the reduced SIR model, using Eq. (11). To be more specific, one can readily solve S(τ peak ) = k = 1 − J decay (τ peak ) with J decay (τ ) from Eq. (71) of [24] for τ peak (for details see Appendix B). In the quoted work, J denotes the cumulative fraction of infected persons, thus J = 1 − S, and J decay applies within the regime of reduced time where the peak time τ peak actually occurs. The resulting time is with in terms of a number of quantities characterizing the solution, that are all expressed in terms of k and η. To be specific, one where W 0 and W −1 are the principal and non-principal solutions of Lambert's equation [23,68]. They are both available like inverse trigonometric functions or elliptic integrals in common software packages (python, Mathematica, matlab, eventually also excel).
For the so-called SK-I approximant, the remaining two parameters b 1 and b 2 left to be specified are given by (Appendix B) while the SK-II approximant is defined by We are not going to interpret all these quantities here, but it may be useful to mention that J ∞ is the finally infected population fraction, and aj max the maximum dimensional rate of newly infected persons (j max is the maximum reduced rate), that occurs at a time that differs from the reduced peak time τ SK peak in Eq. (16) by the two T b terms.

CG approximant
Carvalho and Goncalves [28] (CG) obtained, starting from the reduced SIR model (6) with the help of w = (1 − r)R 0 , for the reduced peak time in I the approximatẽ  14), (c) τ SK−I peak given by Eq. (19), (d) τ SK−II peak given by Eq. (20), and (e) τ CG peak given by Eq. (24). The dimensional peak time t peak is obtained from τ peak upon dividing τ peak by the rate a = β(s 0 + i 0 ). A peak time exists only within the regime η > 1 − k. The corresponding relative deviations between approximant and exact solution are shown in the 2nd row (f)-(j). While the performance of the MT approximant is more accurate than 4.3% over the shown k-η-domain, the SK and CG approximants perform better than MT except within the regime of relatively large η close to 1 − k. This is made more precise in Fig. 2. The color bar for all panels of the first row is shown in (a), the single color bar for the 2nd row is shown in (g). The latter goes from blue (high quality, < 1% deviation from the exact solution) to yellow (low quality, > 5% deviation).
with coefficients given by and where S ∞ is the solution of the nonlinear equation . From Schlickeiser and Kröger [24] we know that the nonlinear equation for S ∞ that remained unsolved in [28] is solved using Lambert's principal function W 0 [23,68] as follows As the reduced timesτ and τ are related by τ = R 0τ according to Eq. (9), with R 0 = a/γ = k −1 according to Eq. (7), the reduced peak time for comparison with the remaining approximants is whereτ CG peak is taken from Eq. (21).

Results and discussion
The exact peak time of the infected compartment calculated from the numerical solution of the SIR equations (4) is shown in Fig. 1a over basically the whole admissible k (horizontal axis) and η (vertical axis) range. Note that we use a semilogarithmic axis and a coloring scheme that reflects log 10 (τ peak ) to appreciate all details, and that the figure is exactly reproduced if we solve Eqs.
(1) or (6) numerically instead. The advantage of (4) or (6) over (1) is, that they have no redundant parameters. The corresponding performance of the four approximants CG, SK-I, SK-II, and MT is shown in Figs. 1b-e, while the relative deviation between approximants and exact solution is given by Fig. 1f- While the performance of the MT approximant is rather accurate over the whole k-η-domain, the SK and CG approximants perform better than MT except within the regime of relatively large η close to 1 − k (in the neighborhood of the upper bound of the colored region, where it transits to white in Fig. 1). This is made more precise later below.
To quantify the performance of the approximants at very small η (at the bottom of Fig. 1d-e), let us mention the following feature, Furthermore, both approximants share the following feature at At the upper bound η = 1 − k the MT approximant correctly vanishes, while the SK approximant does not vanish exactly for all k. This feature causes the SK approximant to become poor in the vicinity of η = 1 − k (see the yellow regime in Fig. 1g-h).
To summarize, Fig. 2 shows the regimes of superior performance of the four approximants. As long as η exceeds a critical η c , the MT approximant is more accurate than the other three. The inequality (27) is thus not fulfilled, and the MT approximant should be used.
To get an impression about the relevant regions in k-η space, Fig. 2 shows in addition many of the published artificial and real data points collected in Tables 1 and 2.

Summary and conclusions
We have compared the quality of the available approximants point-wise, i.e., for each possible case of model parameters separately. Overall, if we limit ourself to the broad regime of k ∈ [0, 1] and η ∈ [10 −9 , 1] the mean relative deviation ⟨∆⟩ between the exact numerical solution and the approximants can be extracted as well. We obtain ⟨∆⟩ ≈ 1.53% for the CG approximant, ⟨∆⟩ ≈ 1.81% for the SK-II approximant, 2.38% for the MT approximant, and a large value of 24% for the SK-I approximant, as it fails dramatically at large η and small k, as can be seen in Fig. 1(g).
Despite this, each of the approximants fails by more than 50% at Fig. 3. Combined approximant, if the best approximant for each region in k-η space is used. This combined approximant has an accuracy smaller than 3% for any choice of SIR parameters. The region of worst performance of the combined approximant is located in the vicinity of k ≈ 0.15 and η ≈ 0.25, while for η < 0.01 or k > 0.5, the performance is just excellent (deviations below 1%). least in some of the k-η region. This fact highlights the necessity to use the best available approximant, given by Fig. 1, depending on the k-η regime that applies to a real situation at hand. Using the combined approximant, the best approximant depending on the (k, η) pair value, the relative deviation to the exact result is quantified in Fig. 3.
There are many examples to which the MT approximant had been applied in the original work [25]. We list them in Table 1. Table 2 confirms that the MT approximant is more accurate when η > η c , while all approximants are good and work very well.
We include in the last three rows of the table a case that is not allowed in the semi-time SIR model, the case of k > 1, as it has been discussed in [25] as well.
Examples of relevance for the ongoing Covid-19 pandemics in 60 countries are available from [17]. The authors analyze the first and second Covid-19 wave in real time. Rather than adding all additional (k, η) pairs from 2021-04-13 to our tables, we have included them all as orange and white circles corresponding to first and second waves to Fig. 2. There is a general trend of a decreasing k = γ /[(s 0 + i 0 )β] for the second pandemic wave, even though the fraction s 0 + i 0 = 1 − r 0 of the remaining unrecovered population has diminished during the second wave. All crosses reside either in the yellow or white regions.
This means, that for most cases of relevance for the ongoing pandemics the existing SK-II and CG approximants are the most accurate. Still, we had to clarify here the relationship between different notation, dimensional versus dimensionless versions of the SIR equations, to allow for a direct comparison. Here, we clarified the correspondence between equivalent versions. The MT approximant might be a good choice if one does not want to calculate a Lambert value, or if the Lambert function is not available at all within the computational environment.
One should keep in mind that the peak time of the daily reported new cases does not coincide with the peak time of the infected compartment. While the former solves d(si)/dt = 0, the latter solves di/dt = 0, and the final fraction of infected population is not related to an integral over i(t), but equals β ∫ ∞ 0 s(t)i(t) dt. To be precise, the peak times differ by the two T b terms in Eq. (16) above. Both terms tend to vanish as k approaches unity. The peak time of the newly infected population fraction is given by τ daily newly = (2/a 3 ) tanh −1 (a 1 /a 3 ) with a 1 and a 3 given by Eq. (18), as shown in [24]. This measurable peak time is thus identical for approximants SK-I and SK-II.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability statement
All data generated or analyzed during this study are included in this published article.

Acknowledgments
We thank the referees for their constructive and helpful comments.
Appendix A. Derivation of τ MT peak (k, η) Our Eq. (14) for the reduced peak time τ MT peak arises from the dimensional peak time denoted by t p2 in Eq. (7b) of [25]. In this latter work one of us (MT) had already shown that his approximant supersedes alternate approximants such as their Eq. (7a) so that there is no need to compare with those. Using the abbreviation S m = γ /β he obtained where we here introduce abbreviations for the otherwise rather lengthy original expression as follows So far we have just reproduced the existing expression for t p2 that had been obtained for the dimensional SIR model (1). To convert this into the dimensionless counterpart, we have to make use of Eq. (8), i.e., we have to replace i 0 , γ , β, and we have to use Eq. (9), which states τ MT peak (k, η) = at p2 . As a result, the s 0 , r 0 , a Table 2 For all entries of Table 1, where the parameters k, η and a of the reduced SIR system can be found, we here compare the performance of the approximants (14), (19), (20), and (24) where we have re-used the quantities A-E defined by (15) to allow for a direct comparison between (A.1) and (14). The s 0 indeed drops out because only the ratio A 3 /A 6 appears in (A.1), and the remaining sign is taken care of the asymmetry of tanh −1 . The rate a still appearing in A 2 drops out in the reduced (dimensionless) peak time because τ MT peak (k, η) = at p2 according to (9). We have thus shown that the reduced time depends only on k and η, and that (14) is the reduced peak time that corresponds to the dimensional peak time t p2 in [25].

Appendix B. Derivation of τ SK peak (k, η)
This appendix provides details on how to read off Eq. (16) from the results obtained by Schlickeiser and Kröger [24]. Equation (71) in [24] for the cumulative fraction J of infected persons after the occurrence of a peak in the differential rate of newly infected persons, within the so-called 'decay' period for reduced times τ ≥ τ * , reads Eq. (54) of [24] (for the special case of c = b, note also that b 3 = |b 1 | as stated after Eq. (61) in [24]) and reproduced here in (17). Furthermore y b = J ∞ as mentioned after Eq. (69) of [24], and τ * is a characteristic reduced time given by Eq. (73) in [24], identical with the first term on the right hand side of (16). The coefficients b 1 and b 2 for the SK-I and SK-II approximants are given by Eqs. (63)-66) in [24] and reproduced in (19) and (20). Analytic results for the characteristic cumulative fraction J 0 , the differential rate of infections at peak time j max , and the cumulative fraction J * and the peak time of j are stated in Eqs. (48), (49), and (62) of [24]. The quantity J decay (τ ) in (B.1) is the cumulative fraction of infected persons at reduced time τ , and thus J decay (τ ) = 1−S(τ ) as S(τ ) denotes the susceptible fraction at reduced time τ . Since we are interested in the peak time of the infected compartment, and because this peak time is delayed with respect to the peak time of the differential rate of newly infected persons, J decay rather than J rise applies; the latter quantity, valid for τ ≤ τ * , had been derived in [24] as well. Making use of (11), one has J decay (τ SK peak ) = 1−k. To be specific, using (B.1), the equation determining τ SK peak becomes }] .

(B.2)
This equation is readily solved for τ SK peak , the result is given by Eq. (16).