Mathematics of Epidemics: On the General Solution of SIRVD, SIRV, SIRD, and SIR Compartment Models

: The susceptible–infected–recovered–vaccinated–deceased (SIRVD) epidemic compartment model extends the SIR model to include the effects of vaccination campaigns and time-dependent fatality rates on epidemic outbreaks. It encompasses the SIR, SIRV, SIRD


Introduction
Compartmental mathematical models are very popular and successful to describe the temporal evolution of pandemic and epidemic outbursts in populations of large size (for reviews see [1][2][3]).Their forecasts on hospitalization and death rates help policy makers to install non-pharmaceutical interventions and/or vaccination campaigns at optimized times.As suitable compartments one introduces the fractions of susceptible persons (S), infected persons (I), recovered persons (R), deceased persons (D), and vaccinated persons (V) that no longer can be infected.Individual time-dependent rates regulate the transition between the different compartments.The temporal evolution of the epidemic is then determined by the ratios k(t) = µ(t)/a(t), b(t) = v(t)/a(t) and q(t) = ψ(t)/a(t) between the recovery (µ(t)), vaccination (v(t)), and fatality (ψ(t)) rates to the infection (a(t)) rate, respectively.By discriminating, e.g., between different age classes in each compartment, these models can be generalized to a much larger number of compartments in order to investigate the effects of pandemic and epidemic outbursts on persons of certain age groups.However, from the health care point of view, to sufficiently provide enough intensive care beds and facilities is independent from the age of the infected persons.
Historically, the first reasonable compartment model was the susceptible-infectedrecovered/removed (SIR) epidemic model [4][5][6][7].This has been refined to include the effect of vaccination campaigns leading to the susceptible-infected-recovered/removedvaccinated (SIRV) epidemic model (see references cited in [8]).The purpose of this manuscript is two-fold.First, we will investigate new classes of exact solutions to the SIRVD and SIRD equations.Secondly, we will also apply the recently developed analytical approximation [8] for the SIR and SIRV models to the more general SIRVD model.These accurate analytical approximations have been derived for all epidemic quantities of interest such as the rate of new infections J(t) and the corresponding cumulative fraction of infections J(t).The main difference between the SIRVD and the SIRV model is the discrimination between recovered and deceased persons by introducing two different compartments.This is necessary as the omicron mutant of COVID-19 has a much smaller (about an order of magnitude) fatality rate than earlier mutants.This gradually changing fatality rate is not accompanied by a corresponding change in the time dependence of the recovery rate.Therefore, a mathematical description with one combined recovery/removed compartment is not sufficiently accurate anymore.
The organization of the manuscript is as follows.In Section 2, we introduce the starting dynamical equations for all considered compartment models both in terms of the real time t with respect to the onset of the pandemic, and the reduced time τ = t 0 dξ a(ξ).It is beneficial for the analysis to express the equations in a form directly involving the observable quantities, such as the cumulative fraction of new infections J(τ) and the cumulative fraction of vaccinated persons V(τ).It is shown that for given reduced time variations of the ratios k(τ), q(τ), and b(τ) the SIRVD and SIRD equations represent complicated integro-differential equations for the rate of new infections j(τ) as well as the cumulative fraction of infections J(τ), or S(τ) and I(τ).However, the integro-differential equations can also be regarded as simpler determining equations for the sum of ratios k(τ) + q(τ) = κ(τ) for given variations of the ratio b(τ) and the fraction S(τ).This new approach is used in Section 3 to derive fully exact analytical solutions for the SIRVD and SIRD models.Especially for the SIRD model, it is an effective new method to construct a special class of exact solutions depending on two parameters which are chosen as the values of the ratio κ(τ) at the start (κ(τ = 0)) and the end (κ(τ = ∞)) of the epidemic outburst.The new method for the SIRD case is illustrated in Section 4 for three different choices of the two parameters including a detailed investigation of the properties of the constructed solutions.
Sections 5 and 6 are concerned with the second main purpose of this manuscript, namely the application of the approximate analytical solution in the limit of small cumulative fractions J ≪ 1 to the SIRVD model.For general reduced time dependencies of the ratios k(τ), q(τ), and b(τ) the time dependence of all quantities of interest is derived in Section 5, whereas in Sections 6.1 and 6.3 two applications are investigated which were inaccessible to analytical treatment before.Of special interest is the calculation of the death rate d(τ) and the corresponding cumulative fraction of deceased persons D(τ).Main differences occur between the considered compartment models, which reflect two alternative points of view.In Section 6.2, we use the analytic solution and its characteristics to obtain all SIRVD model parameters from reported COVID-19 data.
In the SIRVD and SIRD case with a predescribed fatality rate ψ(t), corresponding to the ratio q(τ), the death rate is proportional to the fraction I(t).Moreover, any different reduced time dependencies of the ratios k(τ) and q(τ) correctly enter the dynamical equations.In contrast, in the SIR models no compartment of deceased persons has been considered.Instead, the total fraction of recovered and removed (by death) populations R tot (t) and the summed recovery/removed rate µ tot (t) are used.Then, the solution for the rate of new infections JSIR (t) is employed to calculate an a posteriori death rate as [6] d a−pos = ψ(t) JSIR (t) from a specified fatality rate ψ(t).Of course, this fatality rate can be regarded as part of the summed recovery/removed rate, so that it also enters the dynamical SIR model equations.However, the main difference remains for the calculation of the death rate: in the SIRVD/SIRD cases it is proportional to I(t), whereas in many SIR models it is proportional to J(t).And the temporal dependence of I(t) and J(t) can be different.As we will show, the disparity is most pronounced when the effect of vaccinations is included and/or when the real-time dependence of the fatality rate ψ(t) and the recovery rate µ(t) are different from each other.The a posteriori approach is not necessarily incorrect but has its own justification.It assumes that the probability to die from the virus infection is only determined by being or having been infected with it, and thus is proportional to the rate of new infections J(t).In contrast, in the SIRD and SIRVD models the probability to die is the same on every day being infected and thus depends on the duration of the infection.A summary and conclusion (Section 7) completes the manuscript.
It is appropriate to highlight the important differences to our earlier work [8]hereafter referred to as KS24.KS24 dealt exclusively with the analytical approximation of the SIRV-epidemics model, which has proven the accuracy of the approximate solution by comparing with the numerical solution of the underlying SIRV equations for a few illustrative examples.But no attempt had been made there to practically apply the analytical solution to monitored data from COVID-19 outbursts in order to extract the key parameters of the SIRV model, namely the values of the ratios of the recovery to infection rate and the vaccination to infection rate.Additionally, KS24 does not distinguish between the fraction of deceased and recovered persons, respectively.The additional D-compartment exists only in the SIRVD (and SIRD) models treated in the present work.The dynamical equations for the SIRV and SIRVD models are qualitatively different as the vaccination rate affects the susceptible S-compartment, whereas the fatality rate affects the infected I-compartment, while the summed compartments remain preserved.Ignoring the D-compartment is a significant oversimplification, especially if the effects of vaccinations is taken into account.As demonstrated in the present work, the time dependence of the rates of deceased persons and the rates of newly infected persons is significantly different, which exhibits itself in different values of the respective peak times and ratios of peak intensities.These significant differences, available in analytic form, are useful to apply a novel and powerful diagnostic method to extract the important pandemic parameters of the SIRVD-model.

SIRVD Model
We start with the SIRVD epidemics model described by four transition rates regulating the transitions between the five compartments: the fractions of susceptible persons (S), infected persons (I), recovered persons (R), removed by death persons (D), and vaccinated persons (V) who can no longer be infected.The transition rates in general are time dependent and different from each other.The infection rate a(t) regulates the transition from S → I, the vaccination rate v(t) the transition S → V, the fatality rate ψ(t) the transition I → D, and the recovery rate µ(t) the transition I → R, respectively (Figure 1).The SIRVD equations read: where the dot stands for a derivative with respect to time t.The five fractions obey the sum constraint at all times t ≥ 0 after the start of the wave at time t = 0 with the initial conditions of the so-called semi-time case where η is positive and usually very small, η ≪ 1.

SIRV, SIRD, SIR and SI Models
The SIRVD model includes as special cases the SIRV, SIRD, and SIR epidemics models.For these three simpler models, one introduces the total fraction of recovered and removed (by death) populations and the summed recovery/removed rate (5) by this modification, the five individual compartments of the SIRVD models are reduced to four compartments in the SIRV-description and to three compartments in the SIRdescription, respectively.Consequently, the SIRVD Equations (1a)-( 2) become which, apart from a slight change in notation (R tot and µ tot instead of R and µ before), agree exactly with the earlier considered Equations ( 13)-( 17) in [8].
The SIRD and SIR models ignore the effect of vaccinations so that v(t) = V(t) = 0.The time evolution in the SIRD model is then given by taking the limit v(t) = V(t) = 0 of the Equations ( 1) and ( 2) yielding Likewise, the limit v(t) = V(t) = 0 of the Equation ( 6) provides the SIR equations for the three remaining compartments: For completeness, we discuss the SI model [9] where additionally µ tot (t) = 0 in Appendix A. For µ tot = v(t), the model is known as the SIS model [9].

SIRVD Equations in Terms of the Reduced Time Variable
It is convenient to introduce the reduced time [6] and the dimensionless ratios with this, the SIRVD Equations (1a)-( 2) can be written as By solving Equations (11) numerically for the case of stationary ratios k(τ) = 0.5, b(τ) = 0.01, q(τ) = 0.1, one establishes quantitatively different temporal dependencies for the various fractions of the four different models SIR, SIRD, SIRV, and SIRVD (see Figure 2).Numerical schemes have been developed especially for the SIR model, see the recent works [10,11].We used a variable-step, variable-order (VSVO) Adams-Bashforth-Moulton PECE solver [12] of orders 1 to 13 to produce Figure 2. The highest order used appears to be 12; however, a formula of order 13 is used to form the error estimate and the function does local extrapolation to advance the integration at order 13.This Figure 2 is meant as a motivation for the following analysis, which has the aim to understand the different behaviors of the four models on the basis of analytical calculations.Our analytical study will also be concerned with, in general, time-dependent ratios k(τ), q(τ) and b(τ).Moreover, it will derive new methods to construct special classes of exact analytical solutions.
To this end, it turns out to be important to introduce the rate of new infections, J(t) = a(t)S(t)I(t), which determines in the absence of vaccination the reduction in the susceptible compartment according to Equation (1a).Its dimensionless counterpart, appearing in Equation (11a), is so that J(t) = a(t)j(τ) and j(τ = 0) = η(1 − η).In terms of the rate of new infections j(τ), and the corresponding cumulative fraction of new infections Equations ( 11a) and (11d)-(11f) readily provide at all times Moreover, Equation (11a) yields where we used Equations (11d), (14), and the first Equation (12).Equation ( 14) also provides With Equation ( 15), we obtain for Equations (11c)-( 11e) and implying for the rate of new fatalities The newly infected fraction j(τ), as opposed the currently infected fraction I(τ), is the fraction that is usually reported by health agencies.

Solution of the SIRVD Equations
With the first Equation ( 14) one finds so that Equation (11d) can be written in the form With the initial condition V(τ = 0) = 0 Equation (22) integrates to providing Likewise, Equation (11b) can be written as where we used Equations ( 15) and ( 24).With the initial conditions Equation (25) integrates to a further integration with the initial condition J(τ = 0) = η leads to in terms of S(τ).Furthermore, Equation ( 25) is equivalent to and also with Equation The derived Equations ( 26)-(30) are nonlinear integro-differential equations for the cumulative fraction of new infections J(τ), and S(τ), respectively.The great advantage of the SIRVD equation written in the form (26) is the direct involvement of observable and monitored quantities, such the cumulative fractions of new infections J(t) = J(τ), and of vaccinated persons V(t) = V(τ).We emphasize that Equations (26)-(29) are exact, and hold for all, the SIR, SIRV, SIRD, and SIRVD models.The only difference is that in the SIRDcase J SIRD (τ) = 1 − S SIRD (τ), because here b(τ) = V(τ) = 0, whereas in the SIRVD case J SIRVD (τ) = 1 − S SIRVD (τ) − V SIRVD (τ).If we succeed in either solving the nonlinear integro-differential Equations (26)-(29), or deriving an accurate approximation for their solution, we arrive at a completely analytical solution of all dynamical variables of the SIRVD model and its special cases in dependence on the reduced time τ.
For completeness, we note that inserting the exact solution (23) for V(τ) in Equation (26) leads to the equivalent nonlinear integro-differential Equations for the rate of new infections (31) Before proceeding to the approximation, we note that in terms of the real time, the Equations ( 21)-(26) read, with the help of Equations ( 9), ( 10), ( 12), ( 14), ( 15), (17), (18), and (31), where we introduced the function In some aspects for negligible vaccinations, the resulting SIRD model is simpler than the SIRVD model, e.g., here the simpler relation J(τ) = 1 − S(τ) holds.As a consequence, the general integro-differential Equations (26)-(30) derived before are eased.In Appendix B, we explicitly list the corresponding equations for the SIRD model.
In the following sections, we will derive a class of special exact solutions (Sections 3 and 4) and derive approximate analytical solutions of Equations ( 21) and (25) in the limit of small J(τ) ≪ 1 (Section 5) following our earlier approach [8].This approximation then leads to analytical approximations of all fractions of the SIRVD epidemics model as a function of time for given time dependencies of the three ratios (10).We will prove the accuracy of this approach by comparing it with the exact numerical solutions of these equations for two illustrative examples (Sections 6.1 and 6.3) of the reduced time dependence of these ratios.The proposed analytical approximation is self-regulating as the final analytical expression for the cumulative fraction J ∞ = lim t→∞ J(t) after infinite time allows us to check the validity of the original assumption Before proceeding, we discuss the differences in calculating the fractions of recovered and deceased persons in the different models.

Fraction of Deceased Persons
The main difference between the considered compartment models occurs for the determination of the fraction of deceased persons D(τ) = D(t) and the corresponding death rate d(τ) = dD(τ)/dτ and d(t) = (dD(t)/dt)/a(t) as a function of the reduced and real time.The disparity is most pronounced if the real-time dependence of the fatality rate ψ(t) and the recovery rate µ(t) are different from each other.In the SIRVD and SIRD cases with a predescribed fatality rate ψ(t), corresponding to the ratio q(τ), the death rate is proportional to I(t), as given correctly by Equations ( 18) and (19) and Equations (32e)-(32f).Moreover, the different reduced time dependences of the ratios k(τ) and q(τ) correctly enter the dynamical Equations ( 26)-(29).
However, in the SIR models, no compartment of deceased persons is considered.As mentioned before, see Equations ( 4) and (5), here the total fraction of recovered and removed (by death) populations R tot (t) and the summed recovery/removed rate µ tot (t) are introduced.With the solution for JSIR (t) and J SIR (t) the rate of new fatalities and the fraction of deceased persons are then calculated a posteriori as [6] with α = D ∞ /J ∞ .The difference to Equations (32e)-(32f), which involve and not directly J(t), is obvious.Only if the temporal dependence of I(t) is not drastically different from the temporal dependence of J(t), the a posteriori approach is justified.As has been emphasized before [8], for COVID-19 outbursts in many countries the monitored cumulative fraction of infected persons, J(t) ≪ 1 has been much smaller than unity at all times, so that then in the SIR model indeed I(t) ≃ j(t).
While Equation (34) refers to past work, we will in the following allow for time-dependent fatality rates ψ(t) so that the a posteriori death rate is given by d a−pos (t) = ψ(t) J(t) corresponding to as a function of reduced time.This way, as opposed to the setting (34), the cumulative D a−pos at infinite time must not coincide with D ∞ .Moreover, in the a posteriori approach the influence of the later introduced fatality rate ψ(t) on the dynamical SIR equations is ignored.This is particularly pronounced if the real-time dependence of ψ(t) differs from the real-time dependence of the recovery rate µ(t).We will elaborate on this below with an illustrative example in Section 6.3.

Special Exact Solutions
Equations ( 28)-( 31) are complicated integro-differential equations in the SIRVD and SIRD case, respectively, for given reduced time variations of the ratios k(τ), q(τ) and b(τ).Hovever, they can also be regarded as simpler determining Equations for the ratios k(τ) + q(τ) = κ(τ) for given variations of the ratio b(τ) and the fraction S(τ).This can be used to derive fully exact analytical solutions for the SIRVD and SIRD models.
We know that S(τ) starts from the initial values S(τ = 0) = 1 − η and monotonically decreases to its final non-negative value S(τ = ∞) ≥ 0 after finite or infinite time.We also require, in accord with Equation ( 12) and the initial conditions specified in Section 2.3, that (36) therefore, we adopt the reduced time variation parameterized by S ∞ , τ max , and τ 0 .While lim τ→∞ S(τ) = S ∞ is formally correct, S(τ) may reach zero after finite time τ fin .In that case lim τ→∞ S(τ) = S(τ fin ) = 0, irrespective of the value of the parameter S ∞ , that may therefore take positive or negative values in the expression (37).One has (throughout this work we use the notation tanh −1 (x) = artanh (x)) from Equation (37).Only if the argument of the tanh −1 resides in the interval [−1, 1], a finite τ fin exists.We are going to see that the special choice of S ∞ = 0 has to be treated with care; it will allow us to absorb with Equation (37), along with a particular choice for τ 0 , the analytic solution of the SI model, for which S(τ) reaches S ∞ at τ fin = ∞.Moreover, the initial condition j(0) = η(1 − η) will be used to establish a relationship between S ∞ , τ max , and τ 0 .We adopt without loss of generality positive values of τ 0 > 0. The choice (37) implies With the first Equation ( 15) we obtain for the dynamical SIRVD Equation ( 25) after inserting Equations ( 37) and (39b) this Equation becomes For given reduced time variations b(τ), the combined rate (41) corresponding to the fraction S(τ) in Equation ( 37) can be inferred.In the following, we simplify our analysis to the SIRD model.

SIRD Model
The choice (37) implies for the SIRD model The rate (42a) is of generalized SI-type (see Equation (A5)) with a width τ 0 different from 2 and a general τ max different from ln(1 − η)/η.We thus investigate for which conditions generalized SI-type rates exactly solve the SIRD equations.
. In all four cases, J ∞ ≪ 1, to be more specific: J ∞ = 0.0014 (dotted-dashed), 0.075 (dashed), 0.0028 (solid), and 0.00012 (dotted).As long as S ∞ ̸ = 0, the final value κ ∞ = lim Y→∞ κ(τ) can be also read off immediately from Equation (51).While the last term in Equation ( 51) diminishes with increasing Y due to the leading cosh 2 (Y), the first terms readily evaluate upon replacing tanh(Y) by unity, so that where we made use again of S ∞ from (43).Calculating κ ∞ for the remaining case of S ∞ = 0 requires more care, as the denominator in the last line of Equation (51) does not anymore increase with increasing Y.For S ∞ = 0, one instead finds Note the apparent qualitative difference between the two cases.While

Constraints on the Function κ(τ)
We recall that values of the function κ(τ) smaller than unity describe the epidemic phase where the infection rate a(t) is larger than the sum of the recovery and death rate.Hence, if still enough susceptible persons are available, one expects a rising rate of new infections J(t).In the opposite case of values of the function κ(τ) greater than unity, the sum of recovery and death rates outnumbers the infection rate so that the rate of new infections will decrease in such a phase.Therefore, it is appropriate to start with values of the function κ(τ) smaller than unity at small reduced times.Subsequently, the function κ(τ) approaches the value κ ∞ at infinitely large reduced times, which can be either smaller or greater than unity depending on the chosen values of τ max and τ 0 .However, not all values for τ 0 and τ max are allowed.
The requirement of having a positive κ(0) ≥ 0 leads to this is automatically fulfilled for values of as the right-hand-side of Equation ( 56) is always smaller or equal than unity.For values of τ 0 ≤ 2/(1 − 2η), it is required that The inequality S ∞ ≥ 0, corresponding to J ∞ ≤ 1, is met as long as the following inequalities hold: that follow from Equation (43).This condition (59), the regime between red and blue dashed lines in Figure 4, is not compatible with condition (57) but holds for certain values of τ 0 > 2/(1 − 2η).The condition (58) further guarantees that also S(τ), and not only S ∞ , is non-negative at all times because S ≥ 0 requires, according to Equation (50), which is fulfilled for all times τ because the left-hand side is always smaller than unity while the right-hand side is larger than unity due to condition (58).In cases S ∞ < 0, where the condition (58) is not fulfilled, S(τ fin ) = 0 vanishes at the finite time τ fin already stated in Equation (38).In this case, the chosen solution (37) or (50) can only be used in the finite time interval 0 ≤ τ ≤ τ fin .Equation (51) can also be written as where we used the identity The function (62) is positive for S ∞ ≥ 0 and non-negative for all times in the admissible interval 0 ≤ τ ≤ τ fin .The ratio κ(τ) is negative and contains singularities if it is used in the regime τ > τ fin .

Details of the Construction of the SIRD Solution
From Equation (52) we obtain, for every S ∞ , tanh with positive τ 0 > 0. We note that this Equation has a solution provided because then the right-hand side of Equation (64) lies within the interval [−1, 1].While κ ∞ = 0 for S ∞ = 0, inserting Equation (64) provides for Equation (53), Setting y = 2/τ 0 , Equation (66) leads to the quadratic equation For general and different values of κ(0) und κ ∞ , one can use Equation (67) to determine the implied value of τ 0 = 2/y and with Equation (64) the resulting value of τ max .The knowledge of τ 0 and τ max then determines for the SIRD model the time variation of the ratio κ(τ) for any value of S ∞ and also of the rate of new infections j(τ) in Equation (45).We thus have found an effective new method to construct a special class of exact solutions for the SIRD model which also applies to the SIR model for ψ(t) = 0. We illustrate this new method for several cases of the ratios κ(0) and κ ∞ in the next subsections.
We have thus presented the analytic solutions of the SIRD model for the case of κ ∞ = 0 and arbitrary κ(0).There is one solution for which S approaches zero at infinite time, and there is another solution for which S reaches zero at finite time τ fin (72).In both cases, the characteristic times τ 0 and τ max , as well as S ∞ appearing in the analytic expression (37) for S(τ), are determined by the initial dimensionless rate κ(0) and η = 1 − S(0).

SIRD Solution for Equal
with the two solutions here the function f (K) > 0 is positive for all values of K > 0 and has the slope It therefore decreases for 0 < K ≤ K 0 with K 0 = 1 − (3η/2), where it attains its minimum ).Hence, only the solution is useful as it provides positive values of We first check the requirement (65) reading or The left-hand side of the inequality (82) can be written as and is fulfilled for all values of K ≥ 0. This is obvious for K ≥ K 0 when the right-hand side of Equation ( 83) is non-negative, whereas its left-hand side is negative.For values 0 ≤ K < K 0 Equation (83) reads which is fulfilled.Likewise, the right-hand-side of the inequality (82) can be written as which is fulfilled for values of K ∈ (0, K 0 − η], when the left-hand side of Equation ( 85) is negative while its right-hand side is positive.Recall that the case K = 0 had been treated in Section 4.1.For large values K > K 0 − η squaring the inequality (85) yields We conclude that the constraint (82) can be fulfilled for all values of K ∈ (0, 2K 0 ].In this interval the solution (80) is given by and shown in Figure 6.The first derivative of Equation ( 88) is so that τ 0 (K) attains its maximum value at K = K 0 .Inserting the solution (88) one obtains for Equation ( 64) yielding The solution (92), also shown in Figure 6, is positive for values of 0 < K ≤ 1 − 2η and negative for values of 1 − 2η < K ≤ 2K 0 due to the property tanh −1 (−x) = − tanh −1 (x).
Obviously τ max = 0 for K = 1 − 2η.In order to have positive values of K 0 and the zero 1 − 2η we should consider only values of η < 0.5.
For completeness we note an alternative τ max (τ 0 ).Equating the right hand sides of Equations ( 52) and (53) leads with R ≡ tanh(ρ) = tanh(τ max /τ 0 ) to Upon multiplying with (1 − R)τ 0 one obtains or the quadratic equation which is solved by so that The requirement |R| ≤ 1 is fulfilled for all values of η < 1 if the argument of the square root in Equation ( 96) is non-negative.The last condition leads to which is identical to the maximum value (90) derived before.Proof: given by Equation (37) with S ∞ , τ max , and τ 0 according to Equations ( 43), ( 97) and ( 88).This figure shows τ 0 (80), τ max (92), as well as ρ = τ max /τ 0 (91) versus K.The dotted-dashed vertical line marks K = 2K 0 = 2 − 3η.For K < 1 − 2η, the plus sign applies in Equation ( 97), leading to positive τ 0 and τ max .In that case j(τ) goes through a maximum in the course of positive time τ > 0. For K > K 0 , the minus sign applies in Equation (97).Here, τ max < 0 and j(τ) monotonically decreases at τ ≥ 0. At K = K 0 , τ 0 goes through a maximum versus K, the maximum value τ 0,E ≈ 2/η is given by Equation (90).Within the range of K values where τ 0 and τ max are shown as solid lines, κ(τ) equals K at all times τ ≥ 0 with a precision of less than 0.1%, while within the regions using dashed lines the κ(τ) deviates from K by less than 5%.The analytical solution therefore captures the regime K ∈ [0.8, 2K 0 ] at η = 10 −5 and the regime extends to slightly smaller K for smaller η.
Obviously, positive values of τ fin (K) are only possible for values of K < 1 smaller then unity, as otherwise the argument of the tanh −1 -function is negative.In order for the argument of the tanh −1 -function to be less or equal unity, the condition is required, which is equivalent to Hence, finite values of τ fin (K) only occur for values of K ∈ (0, 0.5].In Figure 8, we show the final time (103) for different values of K below 0.5.The final time is monotonically decreasing from its infinitely large value above K = 0.5 with decreasing values of K.
. While for K < K 0 the κ(τ) varies significantly with reduced time τ, for K > K 0 the form (37) solves the SIRD equations with constant κ(τ) nearly exactly.Lines in (a) are given by Equation (102), which coincides with the numerical solution.Panel (b) offers a zoom into (a) for K 0 < K < 2K 0 .Shown is the numerical solution for κ(τ) − K, which departs from the constant by completely negligible amounts.(c) S(τ) and (d) j(τ).For K < 1/2, S(τ) reaches zero at τ fin given by Equation (110c), the curves for K = 0 (SI model) and K = 1/2 nearly coincide as they have similar τ max and share S ∞ = 0, the curves for K = 1 and K = 1.2 are also indistinguishable by eye and stay at nearly constant S ≈ 1 − η.The bullets mark S(τ fin ) (not S ∞ ) in (c), and j(τ fin ) in (d).
Of particular interest are the cases in Figure 7 where the combined ratio κ(τ) at the start and the end of the epidemic outburst have values below K ≤ 0.5, which corresponds to infection rates being twice as large as the sum of the recovery and fatality rates.In this case, the epidemic outburst, described by the SI-type cosh −2 [(τ − τ max )/τ 0 ]-distribution for the rate of new infections shown in Figure 7d, is so dominated by the rapid infections that at the finite reduced time τ fin the outburst suddenly terminates as with S(τ fin ) = 0 (as seen in Figure 7c) no more persons are available to be infected.The situation is comparable to a smooth-running car where suddenly the car engine stops as no more fuel is available.For values greater than K > 0.5, the epidemic outburst lasts infinitely long.Next, we consider the limits of Equations ( 102) and ( 103) for small values of η ≪ 1.

Limit for Very Small Values of I
and Equation (103).For small values of η ≪ 1 we have to distinguish between the cases (i) where |K − |K 0 | ≤ 2η and (ii) where |K − K 0 | > 2η.We consider both cases in turn.

Case |K − K 0 | ≤ 2η
Here, we approximate for small values of η ≪ 1 Consequently, we obtain whereas τ fin = 0 as all K-values in this interval are greater than 0.5.

Finite Time τ fin (K)
As shown before, c.f. Equation ( 105), finite values of the finite time only occur for values of K ≤ 0.5 < K 0 .Consequently, for small η ≪ 1 we approximate with Equation ( 111) only provides real values if the denominator of the ln-function is positive, leading to the requirement which is consistent with the earlier limit (105) and the time (111) can be further approximated as the function (113) monotonically decreases from its infinite large values above (1 − η)/2 with decreasing values of K in agreement with Figure 8.

Values
With which correctly approaches For values of (1 because for small η ≪ 1 one has K 0 ≃ 1.The solution (118) holds for all times τ in the range (1 For smaller values of K < (1 − 3η)/2 the solution (118) holds only at finite times τ < τ fin (K) given in Equation (113).
As τ 0 has to be positive, the solution (129b) can only be used for values of 1 < κ ∞ < 2 − 3η.

Approximate Analytical Solutions of the SIRVD Model
Recently, it has been demonstrated that an accurate analytical approximation to the SIR and SIRV equations [8] exists if the cumulative fraction of new infections J(t) = J(τ) at all times is much smaller than unity.For the COVID-19 outbreaks in many countries this assumption is well fulfilled.
According to Equation ( 19) in [8] the approximation 1 which approaches where the last approximation holds due to the adopted smallness J ∞ ≪ 1.But we keep the J ∞ in the solution (138) in order not to violate the restriction By integrating Equation (139), the corresponding cumulative fraction is given by Likewise, Equations ( 15) and ( 17)-( 19) provide the approximations and

Difference between the SIRVD Death Rate and the a Posteriori Death Rate
As it marks the difference of the death rate (144) in the SIRVD model to the a posteriori death rate d a−pos (τ) = q(t)j(τ) (35), we consider first the ratio For the SIRD model with negligible vaccinations (b(τ) = 0), the ratio (145) is independent of reduced time and practically equals unity as J ∞ ≪ 1.However, for the SIRVD model including vaccinations the ratio (145) depends on time.The SIRVD death rate has a flatter τ-dependence than the a posteriori death rate.

Conditions for Extrema
As before [8], it is straightforward to calculate the first and second time derivatives of the approximative solution (139) as Consequently, extrema of the rate of new infections occur at reduced times τ E determined by where we indicate that the right-hand side of this Equation is smaller than or equal to unity.Hence, no extrema of infections occur for a sum of variations In the alternative case of reduced time intervals, where extrema are possible.From Equation (146b) one obtains Hence, the extrema are maxima if is positive.They are minima if is negative.There can be multiple minima and maxima depending on the reduced time variation of the ratios k(τ) and b(τ).The extreme values of the rate of new infections are given by All results given in this Section for the SIRVD model include as special cases the corresponding quantities in the SIRD model for b(τ) = 0 and the SIR model for b(τ) = q(τ) = 0, respectively.The excellent agreement with the respective numerical solutions of the SIR and SIRV models [8] has proven the validity of the analytical approximation using the smallness of J(τ) ≪ 1.

SIRVD Model Applications
Next, we illustrate our results from the preceding section with two examples.While Section 6.1 is concerned with the case of stationary ratios, in Section 6.3 we work out the case of a gradually decreasing fatality rate.

Stationary Ratios
We first consider the approximative solutions ( 138) and (139) in the special case of stationary ratios This case corresponds to the numerical solutions of the SIRVD model shown before in Figure 2. If we introduce κ 0 = k 0 + b 0 we may use the earlier results of [8] to obtain and Provided κ 0 + b 0 = k 0 + q 0 + b 0 < 1, the rate of new infections (156) attains its maximum value at the reduced time For the cumulative fraction one finds in terms of the lower incomplete gamma function γ.For infinitely large times, the fraction (159) approaches the final value Likewise, the SIRVD death rate (144) in this case becomes It differs from the a posteriori death rate by the factor e b 0 τ reflecting the ratio between I(τ) and j(τ).The rate (161) peaks at the time which is greater than the peak time τ m of the rate of new infections, since ln(x) is monotonically increasing with increasing x.The maximum SIRVD death rate is smaller than the maximum a posteriori death rate.The corresponding cumulative fraction of fatalities is given by with the final value It has been shown before in Equation (57) of [8] that for small values of b 0 ≪ 1 and b 0 < κ 0 = k 0 + q 0 < 1 the cumulative fractions as well as their final values have to be calculated with the leading order approximation for large values of z ≫ 1: Applying this approximation to Equations ( 160) and ( 165) with z = b −1 0 provides We note that the ratio D ∞ /D a−pos ∞ does not depend on the value of b 0 .
These differences are clearly seen in the panels of Figure 10.The SIR, SIRV and SIRD curves, representing a posteriori death rates, have a significant different reduced time dependence compared to the SIRVD curve.For b 0 ≪ k 0 + q 0 < 1 one infers from Equations ( 159) and (162), For the values of k 0 = 0.5, b 0 = 0.1 and q 0 = 0.1 shown in Figure 10a,b, the maximum death rate is a factor 1.54 smaller than the maximum a posteriori death rate, the final cumulative fraction of fatalities is a factor 1.67 larger than its a posteriori value, and the death rate peaks at the later time τ d = 5.1 as compared to the peak time τ m = 3.6 of the a posteriori death rate.For this example, the exact numerical ratio is τ m /τ d ≈ 0.69, to be compared with ≈0.67 from the analytical expression (168).The performance of the analytical solution j(τ) (156) of the SIRVD model over a wide range of its determining parameters κ 0 and b 0 is analyzed in Figure 11.161) and ( 164), for any k 0 (or q 0 ). (a) shows the range κ 0 , b 0 ∈ [0, 2], while (b) offers a zoom into the bottom left region κ 0 ∈ [0, 0.5] and b 0 ∈ [0, 0.2] of pronounced relative deviations, while the absolute deviation is still small, as visible from Figure 10e,f, where we show a case that resides in the yellow region.

Application to Measured COVID-19 Data
Here, we exemplify how to obtain all SIRVD model parameters from measured COVID-19 data of a completed pandemic wave if all rates are considered time independent.Using results from the foregoing section, we are going to demonstrate how to extract all parameters analytically from a few characteristics contained in the reported data.
To this end, we collected data capturing the first pandemic COVID-19 wave for several countries.Usually, only the cumulative number of newly infected persons, N J(t), and the cumulative number of deceased persons, ND(t), are reported, while the population size N can be considered known.Sometimes, R(t) is reported as well, but is typically useless, as not all recovered persons report their recovery, and R(t) is then simply estimated based on J(t) and D(t).Note that the current fraction of infected population, I(t), as well as the susceptible population fraction, S(t), are usually not measurable.
We thus rely on the reported J(t) and D(t) time series as a function of time t, as well as their derivatives J = dJ/dt and D = dD/dt, starting at the time t = 0 of the outbreak, which we define here to coincide with the day at which 10 persons have been infected so far, N J(0) = 10.The most robust quantities that can be obtained from such time series are the final plateau values J ∞ and D ∞ , when J and D approached zero or a value that is very small compared with the peak values Jmax and Dmax of the measured rates.As proven in Appendix D, the following procedure can be followed: (i) calculate the ratio G ∈ (0, 1) defined in Equation (A41), i.e., there is no set of positive rates b 0 and κ 0 that allow one to capture the data subject to the inequality b 0 + κ 0 < 1, as explained in Appendix D. If, however, G ≤ e −1 , calculates a real-valued positive U > 1 with the help of the principal branch W 0 of Lambert's W function, where the choice of principal versus non-principal branch depends on the magnitude of G.While Lambert's W function is routinely available in engineering software, Figure 12 can alternatively be used to look up U for given G, (iii) with U and the three measured ratios t d /t m , D ∞ /J ∞ , Dmax / Jmax at hand, calculate according to Equation (A41); note that U > 1 implies κ 0 < 1, (iv) obtain the three dimensionless SIRVD rates via Equation (A37), i.e., Finally, the dimensional rates, if needed, are given by Equations ( 10) and (158), since a 0 = τ m /t m .This procedure has been followed to create Table 1 for a few selected countries, for which the assumptions b 0 > 0, κ 0 > 0, and κ 0 + b 0 < 1 underlying the present treatment apply, using public available data [13].With the parameters at hand, the time-evolution not only of the reported, but also of the remaining compartments S, R, and V is predicted by the SIRVD model.The outlined simple recipe must not produce the best possible fit to the measured data by any measure, but it ensures that the three characteristic values are exactly reproduced.Moreover, in an analogous fashion it is possible to use the analytic solution to efficiently extract the SIRVD parameters from measured data well before peak times have been reached, for example, using polynomial coefficients of the time series, ratios between cumulative fractions at different times, etc.This allows to forecast the evolution of all compartments.

Gradually Decreasing Fatality Rate
Here, we investigate the approximative solutions (138) and (139) in the special case of stationary ratios k(τ) = k 0 , b(τ) = b 0 but the decreasing fatality rate with constants q 0 , q 1 and G.The fatality rate (174) starts from the constant value q 0 and decreases with the typical time scale G −1 to its final constant value q 1 .Such a behavior accounts well for the COVID-19 omicron which had a much smaller (about an order of magnitude) fatality rate than the earlier mutants occurring about a year earlier.Such a slow gradual decrease is well captured by the adopted fatality rate (174) in the case G ≪ b 0 < 1, allowing us the linear approximation if necessary.The vaccination rate (155) is independent from the fatality rate q(τ) and also applies here, whereas with the rate (174) the rate of new infections (139) becomes where we made use of the abbreviation Consequently, one finds for Equations ( 141)-( 144) and respectively, whereas the cumulative fraction of infections is given by The analytical expressions are in excellent agreement with the numerical solutions of the SIRVD model in the presence of time-dependent q(τ).In Figure 13, we show selected results for j(τ) and d(τ), while all other quantities are equally well captured.
Both Equations ( 181)-( 182) are nonlinear and cannot be solved in closed form (this is known from the works of Évariste Galois who died during a duel at the age of 20 in 1832.His important contribution, known as Galois theory, was recognized and published 14 years later by Joseph Liouville [14]).However, if we use, as in Equation (175), the approximation e −Gτ m.d ≃ 1 − Gτ m,d Equations ( 181)-(182) can be approximated as and This latter Equation (184), treated below, is of the form of Wright's transcendental equation.It is easy to see that Equations (176)-(184) in the limit q 1 = 0 correctly reduce to the earlier results in Section 6.1.The applicability of Equations ( 183) and (184) requires G ≪ b 0 .Introducing the positive z m = b 0 τ m and the positive combination a m = b 0 C 0 /q 1 G of parameters, the first Equation ( 183) reads e −z m = C 0 [1 − (z m /a m )] and thus determines z m in terms of a m and C 0 .This latter equation is solved in terms of the non-principal Lambert function W −1 (see Appendix G in [15]) as where we have just re-inserted z m and a m in the second part of Equation (185).
As detailed in Appendix C, the solution of Equation ( 184) is more involved.As for the calculation of τ m , it is helpful to introduce a dimensionless time z d = b 0 τ d , and to simplify Equation (184) using an appropriate combination of the five semipositive parameters k 0 , b 0 , q 0 , q 1 , G, and z d .As shown in Appendix C, it is possible to write Equation (184) in the form of Wright's transcendental equation, or equivalently a = −X − b coth(X/2) (187) upon introducing X, a, and b via Since q 1 < q 0 , to ensure positive q(τ), one has b > 0. Further, a < 0 since G ≪ b 0 .
For positive values of b and negative values of a, no solutions to Equation (187) with negative X are possible.This is most easily seen by rewriting this equation as a + X = b coth(X/2).For negative X, the right-hand side is negative, while the left-hand side is positive.Equation (187) hence implies X > 0, and z d > 0 implies X > ln C 1 .Equation ( 186) is equivalent to Wright's transcendental equation and known [16] to exhibit real-valued solutions X only for a limited range of a and b values.Using the x-parametric form of the envelope a = −x − sinh(x) and b = −1 + cosh(x) [16] we derive the condition for the existence of a τ d value.The a max monotonically decreases with increasing b, starting from a max = 0 and b = 0.
In Appendix C, we obtain the following approximations for the exact solution of Equation (187): Note that X 0 approaches unity for a = a max and b → ∞, while X 2 → ∞ behaves correctly in this limit.As detailed in Appendix C, X 0 ist most precise for sufficiently small a well below a max , X 1 performs extremely well (Figure 14) over the whole a-b range except very close to a = a max and large b, while X 2 has advantages only in the regime where X 1 fails.To summarize, τ d exists as long as a < a max , and is then well approximated by with a and b expressed in terms of k 0 , b 0 , q 0 , q 1 , and G via Equation (188).In Figure 13, we mark the analytical τ m and τ d for a few cases.As visible, they capture the peak times of the analytical solutions for j(τ) (176) and d(τ) (178c), which moreover are indistinguishable from the numerical solution.Since it is impossible to visualize the performance of the numerical versus analytical results for the case of time-dependent q(τ), characterized by 6 parameters k 0 , b 0 , q 0 , q 1 , G, and η, we have randomly chosen sets of parameters with k 0 , q 0 , b 0 ∈ [0, 1], q 0 ∈ [0, Since it is impossible to visualize the performance of the numerical versus analytical results for the case of time-dependent q(τ), characterized by 6 parameters k 0 , b 0 , q 0 , q 1 , G, and η, we have randomly chosen sets of parameters with k The number of sets is (a) 100 and (b) 5000.In (a) we visualize the performance of analytical versus numerical results for τ m (185) and τ d (191), while (b) shows τ d versus τ m , both numerical valus (red circles) and analytical expressions (black squares).The numerical results suggest τ d ≥ τ m , while this is not reflected by 1.5% of the analytical results.

Maximum Rate of New Infections and Death Rate
With the peak times (185) and (191) inserted in Equations ( 176) and (178c), respectively, the maximum rates of the new infections and the maximum death rates are given by where we made use of the determining Equation (181), and
We observe earlier peak times for the a posteriori death rates, following the rate j(τ) of newly infected individuals, compared with the SIRVD death rate d(τ), following I(τ).Additionally, the SIRVD death rates are significantly larger compared with those obtained through the a posteriori treatment, with direct implications for the corresponding D ∞ values.Our predictions for these differences are suitable for being corroborated with past monitored pandemic data on the rate of new infections and the death rates of sufficient quality with negligible underestimation, i.e., negligible dark numbers.

Summary and Conclusions
We have investigated a special class of exact solutions as well as accurate analytical approximations of the SIRVD and SIRD compartment models.For nonlinear models with a high-dimensional parameter, space analytical solutions, exact or accurately approximative, are of high importance and interest: not only as suitable benchmarks for numerical codes, but especially as they allow us to understand the critical behavior of epidemic outbursts as well as the decisive role of certain parameters [17].For given reduced time variations of the ratios k(τ), q(τ), and b(τ), the SIRVD and SIRD equations represent complicated integrodifferential equations for the rate of new infections j(τ) as well as the cumulative fraction of infections J(τ), or S(τ) and I(τ).However, the integro-differential equations can also be regarded as simpler determining equations for the sum of ratios k(τ) + q(τ) = κ(τ) for given variations of the ratio b(τ) and the fraction S(τ).This new approach has been used to derive fully exact analytical solutions for the SIRVD and SIRD models.Especially for the SIRD model it is an effective new method to construct a special class of exact solutions depending on two parameters which are chosen as the values of the ratio κ(τ) at the start (κ(τ = 0)) and the end (κ(τ = ∞) of the epidemic outburst.The new method for the SIRD case is illustrated in Section 4 for three different choices of the two parameters including a detailed investigation of the properties of the constructed solutions.Of particular interest are the cases where the combined ratio κ(τ) at the start and the end of the epidemic outburst have the same value below K ≤ 0.5 which corresponds to infection rates being twice as large as the sum of the recovery and fatality rate.In this case, the epidemic outburst, described by the SI-type cosh −2 [(τ − τ max )/τ 0 ]-distribution for the rate of new infections, is so dominated by the rapid infections that at the finite reduced time τ fin the outburst suddenly terminates, as with S(τ fin ) = 0 no more persons are available to be infected.The situation is comparable to a smooth-running car where suddenly the car engine stops as no more fuel is available.For values greater than K > 0.5, the epidemic outburst lasts until infinitely large times.
In the second part of the manuscript, the recently developed analytical approximation [8] for the SIR and SIRV models are applied to the more general SIRVD model.In the limit of small cumulative fractions J ≪ 1, which very often is fulfilled, this approximation provides accurate analytical expressions for all epidemic quantities of interest such as the rate of new infections J(t) and the fraction I(t) of infected persons.One of the referees has kindly informed us that our analytical approach is related to the Perov's fixed point theorem [18].As an aside, when determining τ d in Appendix C we provided accurate approximative solutions to Wright's transcendental Equation (A16), or equivalently, Equation (186).The main difference of the SIRVD to the SIRV model is the discrimination between recovered and deceased persons by introducing two different compartments which affects the calculation of the death rate.In the SIRVD/SIRD cases it is proportional to I(t), whereas in many SIR models in an a posteriori approach it is proportional to J(t).It is shown that the temporal dependence of I(t) and J(t) are different when the effect of vaccinations is included and/or when the real-time dependence of the fatality rate ψ(t) and the recovery rate µ(t) are different from each other.We illustrate these pronounced differences with two applications: one for stationary ratios k 0 , b 0 and q 0 , and one for stationary ratios k 0 , b 0 but a gradually decreasing fatality rate.In our analysis, we observe earlier peak-times for the rate of newly infected individuals compared with death rates.Additionally, our findings indicate significantly smaller death rates compared with those obtained through the a posteriori treatment.We are hopeful that these differences can be corroborated with past monitored pandemic data on the rate of new infections and the death rates of sufficient quality, with negligible underestimation.We did not account for the widely acknowledged 7-day delay between the onset of infection and the resulting death.This delay impacts both the death rate in our SIRVD model and the subsequent a posteriori death rate analysis in a consistent manner.As a result, the derived difference in peak times remains unchanged.The case of stationary ratios allows one to construct a new powerful diagnostics method to extract analytically all SIRVD model parameters from measured COVID-19 data of a completed pandemic wave.The new diagnostics method is applied to the monitored COVID-19 data in four countries.Potential future work on this subject should include (1) incorporating in the SIRVD model spatially heterogeneous situations by adding spatial diffusion, (2) a detailed testing of the predictions with suitable data from past COVID-19 waves also for time-dependent ratios k(τ), q(τ), and b(τ), and (3) the derivations of accurate mathematical approximation for more complicated time variations of the ratios k(τ), q(τ) and b(τ).
Here, we depart from Wright's arguments and solve Equation (A33) with asymptotic expansions of the coth-function.For values of |Z| ≤ 2 we approximate coth(Z/2) ≃ 2/Z[1 + (Z 2 /12)], so that in this range Equation (A33) reduces to the quadratic with the two solutions The term under the square root is not strictly semipositive for any a < a max , but Z ± is real-valued for a < −2 b(6 + b)/3, which is only slightly below a max ; for b < 10 it is only 3% below a max , giving rise to the tiny white region in top right corner of   189), approximants X 0,a , X 0,b , X 1 , and X 2 given by Equation (190).For any X and b, the exact solution corresponds to a = −X − b coth(X/2) according to Equation (187).Two contourlines are highlighted in each panel, if they exist: X = 1 (dashed black) and X = 2 (solid black).In the X 1 panel, the X 1 is undefined in the small white region in the top right corner, defined by a < −2 b(6 + b)/3.In the X 2 panel, the X 2 cannot be evaluated numerically in the yellow region which is characterized by X > 1 and a > a max /2.Bottom row from left to right: The relative deviation ∆X between approximants X 0,a , X 0,b , X 1 and X 2 and the exact numerical solution X, again versus a max /a and b/(b + 1).The color scale ranges from blue (relative deviation below 1%, to yellow (relative deviation above 10%).Approximant X 1 performs very well except in the tiny top right corner (large b, a close to a max ), where X 2 has some advantage.

Appendix D. Recipe for Stationary Rates
In Section 6.2, we showed how to extract all SIRVD model parameters from reported data.Here, we derive the required relationships.For stationary infection rates a(t) = a 0 the reduced time is defined by τ = a 0 t as we have adopted without loss of generality t = 0

Figure 11 .
Figure 11.Performance of the analytic solution of the SIRVD model for η = 10 −5 and constant rates in κ 0 -b 0 space.Shown is the mean percentage of deviation between the analytic j(τ) from Equation (156), and the numerical solution.The plot looks nearly identical for all remaining differential and cumulative quantities, including d(τ) and D(τ) from Equations (161) and (164), for any k 0 (or q 0 ). (a) shows the range κ 0 , b 0 ∈ [0, 2], while (b) offers a zoom into the bottom left region

Figure 12 .
Figure 12.The nontrivial solution U(G) ̸ = 1 solving U ln(UG) = ln G for G ∈ (0, 1) can be expressed in terms of Lambert's W function, c.f., Equation (170).Depending on the value for G, the principal branch W 0 (black) or non-principle branch W −1 (red) applies, as shown.The black-red crossover is at G = e −1 ≈ 0.3679 and U = 1.The function U(G) is used to map characteristics of the measured data to the SIRVD model parameters k 0 , q 0 , b 0 , and a 0 , or also its dimensional counterparts µ 0 , v 0 , and ψ 0 , as described in Section 6.2.(a) linear-linear and (b) double-logarithmic show the same function U(G).

15 Figure 14 .Figure 15 .
Figure 14.(a) The exact solution X to Equation (186) or (187) versus a max /a and b/(b + 1) with the negative a max defined by Equation (189).(b) Approximant X 1 given by Equation (190b).It is undefined in the small white region in the top right corner, characterized by a < −2 b(6 + b)/3.Two contourlines are highlighted in each panel: X = 1 (dashed black) and X = 2 (solid black).For the remaining approximants X 0 and X 2 see FigureA1in Appendix C.

Figure A1 .
For this reason the question about the suitable sign in Equation (A35) cannot be answered by considering the limit b → 0 and a → a max , as before.Instead, consider the following: Inserting a = Z + χ sinh(Z) and b = χ[cosh(Z) − 1] solves Equation (A33) for any choice of χ , i.e., one can writea = Z − b(b + 2χ), Z = − cosh −1 b + χ χ , (A36) for positive χ.Using χ = 1 we had derived a max .Using χ = 0 in Equation (A36) yields b = 0 and a = Z, while the second relation in Equation (A36) in this limit is just the identity Z = Z.Inserting this special case of a = Z and b = 0 into Equation (A35) yields Z ± = (Z ± Z)/2.As we expect to find Z ± = Z, the suitable solution is Z + .With X 1 = −Z + < 2 this completes the derivation of X 1 in Equation (190b).

Figure A1 .
Figure A1.Top row from left to right: The exact solution X to Equation (186) or (187) versus a max /a and b/(b + 1) with the negative a max defined by Equation (189), approximants X 0,a , X 0,b , X 1 , and X 2 given by Equation (190).For any X and b, the exact solution corresponds to a = −X − b coth(X/2) according to Equation (187).Two contourlines are highlighted in each panel, if they exist: X = 1 (dashed black) and X = 2 (solid black).In the X 1 panel, the X 1 is undefined in the small white region in the top right corner, defined by a < −2 b(6 + b)/3.In the X 2 panel, the X 2 cannot be evaluated numerically in the yellow region which is characterized by X > 1 and a > a max /2.Bottom row from left to right: The relative deviation ∆X between approximants X 0,a , X 0,b , X 1 and X 2 and the exact numerical solution X, again versus a max /a and b/(b + 1).The color scale ranges from blue (relative deviation below 1%, to yellow (relative deviation above 10%).Approximant X 1 performs very well except in the tiny top right corner (large b, a close to a max ), where X 2 has some advantage.

Table 1 .
[13]r part of the table: Quantities extracted from reported COVID-19 data for various countries (population size N).Time of outbreak defining t = 0 (day in 2020), cumulative infected persons Nη at outbreak, peak times t m (infected) and t d (deceased), corresponding peak amplitudes Jmax and Dmax as well as final fractions J ∞ and D ∞ obtained from reported time series J(t) and D(t) of the cumulative fraction of infected and deceased population.The differential counterparts, dimensional rates J(t) = dJ(t)/dt = a 0 j(τ) and D(t) = dD(t)/dt = a 0 d(τ) with τ = a 0 t can be extracted directly from the data that is usually provided daily (dt equals one day).The lower part of the table displays dimensionless ratios extracted from the upper part, the resulting three dimensionless rates k 0 , q 0 , b 0 along with the dimensional rate a 0 of the SIRVD model.Results are quite insensitive to the precise choice of the outbreak time.Additional COVID-19 data available from[13].