Analytical Modeling of the Temporal Evolution of Epidemics Outbreaks Accounting for Vaccinations

With the vaccination against Covid-19 now available, how vaccination campaigns influence the mathematical modeling of epidemics is quantitatively explored. In this paper, the standard susceptible-infectious-recovered/removed (SIR) epidemic model is extended to a fourth compartment, V, of vaccinated persons. This extension involves the time t-dependent effective vaccination rate, v(t), that regulates the relationship between susceptible and vaccinated persons. The rate v(t) competes with the usual infection, a(t), and recovery, μ(t), rates in determining the time evolution of epidemics. The occurrence of a pandemic outburst with rising rates of new infections requires k+b<1−2η, where k=μ(0)/a(0) and b=v(0)/a(0) denote the initial values for the ratios of the three rates, respectively, and η≪1 is the initial fraction of infected persons. Exact analytical inverse solutions t(Q) for all relevant quantities Q=[S,I,R,V] of the resulting SIRV model in terms of Lambert functions are derived for the semi-time case with time-independent ratios k and b between the recovery and vaccination rates to the infection rate, respectively. These inverse solutions can be approximated with high accuracy, yielding the explicit time-dependences Q(t) by inverting the Lambert functions. The values of the three parameters k, b and η completely determine the reduced time evolution of the SIRV-quantities Q(τ). The influence of vaccinations on the total cumulative number and the maximum rate of new infections in different countries is calculated by comparing with monitored real time Covid-19 data. The reduction in the final cumulative fraction of infected persons and in the maximum daily rate of new infections is quantitatively determined by using the actual pandemic parameters in different countries. Moreover, a new criterion is developed that decides on the occurrence of future Covid-19 waves in these countries. Apart from in Israel, this can happen in all countries considered.


Introduction
In December 2020 the effective mRNA-based Covid-19 vaccine by the companies Pfizer-BioNTech and Moderna became available. This has lead to intensive vaccination campaigns in many countries over the world with different speeds. As two shots per person are needed for nearly 100 percent protection as of February 10, 2021, Israel with a time t-dependent vaccination rate of v(t) 7.0 × 10 −3 day −1 has the highest, followed by the United Arab Emirates v(t) 4.7 × 10 −3 day −1 , United Kingdom with v(t) 2.0 × 10 −3 day −1 , whereas Germany's vaccination rate v(t) 4.2 × 10 −4 day −1 is significantly smaller.
It is the purpose of the present manuscript to investigate analytically and quantitatively for a given ratio b(t) = v(t)/a(t) of the vaccination rate to infection rate a(t) the effect on the time evolution of the ongoing epidemic waves. We base the analysis on the susceptible-infectious-recovered/removed (SIR) epidemic model [1] augmented by the appropriate vaccination rates leading to the susceptible-infectious-recovered/removedvaccinated (SIRV) epidemic model. This model is a dynamical system for time-dependent quantities S(t), I(t), R(t) and V(t) denoting the relative fractions of currently susceptible (S), infectious (I), recovered/removed (R) and vaccinated (V) persons in the considered population of N persons as a function of time t (Fig. 1). In the case of negligible vaccination, assuming a constant ratio between infection and recovery rate, considerable improvements on the analytical modelling of epidemics with this compartment model has been achieved citeks20,sk21. We make frequent use of these improvements in the following. Application of this earlier work [2] to the monitored second waves in many countries has indicated initial infection rates of the order of a(0) ∈ [0.1-1.0] day −1 , considerably greater than the above vaccination rates. However, it is important to emphasize an essential difference: whereas the initial vaccination and the recovery rate can be directly used to estimate the corresponding typical time scales T v 1/v(0) and T r 1/µ(0) for vaccinations and recovery, respectively, the initial infection time scale T i = 1/[a(0)I(0)] additionally depends on the initial fraction I(0) of infected persons at the onset of the 2nd wave. With this, and depending on the country, the vaccination time scale T v is often comparable with the infection time scale T i .
The inferred infection rates are slightly larger than the initial recovery rates µ(0) so that the ratio of the two k = µ(0)/a(0) ∈ [0. 8,1). This is consistent with the result [3] that for a pandemic outburst with a prominent peak at a later time the ratio k has to be smaller than k < 1 − 2η, where η = I(0) is used. For most second waves η 1 is negligibly small, so that the determined values of k less than unity are fully consistent. Then with the above noted vaccination rates the ratio of the vaccination to initial recovery rates b = v(t)/µ(0) is considerably smaller than the ratio k: we expect values of b ∈ [5 × 10 −4 , 10 −2 ] k. Consequently, the difference α = k − b is positive and only slightly smaller than k. We may refer to α as the effective ratio as compared to the ratio k.
It is the purpose of this manuscript to determine quantitatively the influence of this small reduction of the ratio k on the time evolution of the pandemic wave. As in earlier work before our analytical calculations are based on the assumption that the ratios k and b are constants; but they hold for arbitrary time dependent infection rates a(t), where, however, the recovery and vaccination rate follow the very same time dependence as a(t).
Compartment models for epidemics with vaccinations have been considered before [4][5][6][7] but they were more concerned with optimizing the control of the epidemics by vaccination with limited resources. Our main goal of this study is to derive analytical solutions for the dynamical SIRV to be presented in Eqs. (1)- (5). In order to keep the analysis as simple and transparent as possible we ignore complicating issues such as age grouping, vital dynamics and/or spatial spread effects that have been recently investigated in the literature [8,9] with numerical solutions.

General SIRV equations
The SIRV-model is a dynamical system for time-dependent population fractions S(t), I(t), R(t) and V(t) introduced above. The fractions add up to unity, S(t) + I(t) + R(t) + V(t) = 1 at all times (sum constraint). The SIRV differential equations accounting for vaccinations of the susceptible persons with the vaccination rate v(t) reaḋ where the dot denotes a derivative with respect to time t. The SIRV equations are supplemented by semi-time initial conditions (defining t = 0) with η ∈ (0, 1) denoting the fraction of infected persons at time t = 0. One of the four dynamical equations can also be replaced by the sum constraint. These initial conditions are sufficient to capture applications of the SIRV equations with non-vanishing R(0) and V(0) as this is identical with the present initial condition upon subtracting the nonvanishing ∆N = N[R(0) + V(0)] from N. The resulting SIRV quantities are then fractions of the reduced, susceptible or currently infected population, and the fractions S and I with respect to the total population N are obtained by multiplying them both with 1 − (∆N/N). Similarly, the total number of recovered and vaccinated persons is (N − ∆N)R plus NR(0), same for V. More formally, ifX denotes one of the SIRV quantities that fulfills initial conditionX(0), the time-evolution ofX is given by the time-evolution of X viaX = χX for X ∈ {S, I} andX =X(0) + χX for X ∈ {R, V} with χ =Ĩ(0) +S(0).
The initial conditions (5) are specified byη =Ĩ(0)/χ, and the parameters to be introduced next have to be adjusted as well,k = χk andb = χb.
The fractions S(t), I(t), and R(t) are usually not measurable with high confidence, while the daily new number of infected persons, denoted bẏ J = a(t)SI, (6) and the fraction of vaccinated persons, V(t), are two quantities that can be more easily measured, and are usually reported. Using Eqs. (1) and (4) the total cumulative number fraction J of persons, that have been infected up to the time t, is related to the SIRV quantities via or equivalently, upon making use of the sum constraint, by the sum of currently infected and currently recovered, Following previous works, and to make sure that the SIRV model has any predictive power, we will allow for arbitrary time-dependent rate a(t) throughout this work, but assume at the same time that the remaining rates µ(t) and v(t) exhibit the identical dependency with respect to time. This leaves us with two time-independent model parameters, As we will demonstrate the two parameters k and b together with the initial fraction of infected persons η completely determine the temporal evolution of the pandemic wave in the reduced time τ = t 0 dξa(ξ).

Condition for pandemic outburst
It is instructive to calculate with the first two SIRV equations the initial variation of the daily number of the newly infected population fraction, where we have omitted the argument t for all functions, and used Eq. (8). In order for a pandemic outburst with initially rising rates of new infectionsJ(0) > 0 to occur the bracket on the right-hand side of Eq. (9) has to be positive at the starting time t = 0, so that with initially constant rate valuesȧ(0 Making use of the first two initial conditions (5) the outburst condition reads As k and b have positive values the condition (10) implies: (i) If initially more than 50 percent (η > 0.5) are infectious, no new pandemic outburst will occur. However, such high values of η are unlikely and unrealistic. (ii) For small given values of η 0.5 and the ratio of recovered to infection rate k, new emerging outbreaks can be fully prevented for values of the ratio of vaccination to infection rate b > 1 − k − 2η 1 − k. The more pathogenic a virus mutation is, the smaller and closer to zero is the value of the ratio k so that the lower limit for b has to be close to unity to prevent a new outburst. (iii) For any finite value of η for modeling epidemic outbreaks the relevant range of the two parameters k and b is 0 ≤ b + k < 1. (iv) As an aside we note that Eq. (10) demonstrates that in the SIR-model with b = 0 no pandemic wave can occur if the parameter k = 1 equals unity. The SIR-model correctly indicates that epidemic waves end in the case k = 1. Therefore the recent criticism [10] on the SIR-model is inappropriate and misguided.

Reduced time
We consider the case where the ratios of the recovery to infection rate, µ(t)/a(t) = k, and the vaccination to infection rate, v(t)/a(t) = b, are semipositive constants independent of time. This assumption still allows us to account for any given time-dependence of the infection rate with the caveat that the recovery and vaccination rate have exactly the same time dependence as the infection rate apart from their different initial values. The introduction of the reduced time scale reduces the SIRV equations (1) -(4) to Eqs. (12) and (13) providing for Eqs. (14) and (15) Both equations integrate immediately to where the integration constants have been determined with the initial conditions (5) holding now at τ = 0. In terms of the reduced time τ the differential new number of infected while the corresponding cumulative fraction is At τ = 0, the cumulative J(0) = I(0) + R(0) = η. In reduced time Eq. (9) corresponds to We have so far expressed R, V, J, and j in terms of S and I. An additional relationship between S and I is provided by the sum constraint, as shown next, and this will allow us to come up with a closed equation for a single variable, to be developed in the next section.
Inserting Eqs. (16)-(17), (20) and (21) provides for the sum constraint In the present manuscript we will derive analytical solutions of Eq. (25) for general non-zero and different values of b = k. We will obtain an implicit analytic solution that expresses the reduced τ in terms of a parameter ψ, while all SIRV-functions including j are expressed in this parameter. We further present a highly accurate analytical approximation for all SIRV-functions as a function of τ. Our new analytical solutions reduce in the appropriate limit to the earlier [2,3,11] solutions for the non-vaccination case b = 0. We also consider as special cases the non-recovery case k = 0 and the special case of equal values of k = b.

Summary of results
Because the following derivation of the solution of the SIRV equations (11)-(15) in reduced time with parameters k and b, and subject to initial conditions I(0) = 1 − S(0) = η and R(0) = V(0) is rather lengthy, we begin by stating the final result. We are able to derive an implicit exact solution τ = τ(ψ) parameterized by ψ, while all SIRV quantities can be expressed in ψ as well. Provided the reduced vaccination rate b exceeds a critical b c , for which we provide the explicit expression (91), we show that the explicit solution of the SIRV equation can be written as follows. With the help of ψ(τ) = ln[S(τ)/I(τ)], we obtain where the differential j and cumulative fractions J of infected persons are given by If b is smaller than the critical b c , provided by Eq. (91), the SIRV model approaches the SIR model that has been treated before, and the SIRV quantities are well captured by a simple linear superposition of the SIR result with the above SIRV solution evaluated at b = b c . Next we derive these expressions and provide additional features of the solution such as the final values at t → ∞, treat special cases such as k = b, and discuss the features of the SIRV equations that gave rise to Eq. (26). The solution holds for any semipositive values of k and b.

Two useful functions
We start the mathematical analysis by introducing the function implying with the usually positive (for initial fractions η < 0.5 of infected persons) initial value The first derivative of the ψ (31) is given by where we used Eqs. (12) and (13) and where the following abbreviation is introduced Consequently, the function Similarly, Eq. (20) is identical with S(τ) = (1 − η) exp[−bτ − R(τ)/k], so that with b > 0 and R residing in the finite interval [0, 1], one has S ∞ = S(τ = ∞) = 0. Consequently, Φ ∞ = 0 ultimately vanishes as well. In this limit Eq. (36) says that dψ/dτ is reaching the constant α, and thus ψ ∞ = ψ(τ = ∞) = sign(α)∞ is infinitely large, while its sign is given by the sign of α, as long as α = 0. For α = 0 the ψ ∞ reaches a constant value.
The initial slope of dψ/dτ evaluates to and is thus negative for all α = k − b < 1 and positive for α > 1. This implies, that the function ψ decreases initially, undergoes a minimum ψ m at Φ m = α before increasing to its final value ψ ∞ = ∞, when α ∈ (0, 1). For negative α, ψ monotonously decreases, while for α > 1 the ψ monotonously increases. For all non-zero values of k and b the function Φ = I + S, however, decreases at all times from its initial maximum value Φ(0) = 1 as its derivative is always seminegative, i.e.

Mathematical analysis
With Eq. (32) the two Eqs. (36) and (38) yield The combination of these two Eqs. (39) then provides or equivalently, with the function which we prefer to write as the derivative of B given by With the function (43) we obtain for Eq. (41) multiplied with (dψ/dτ) = α − Φ from Eq. (36) with the first integral where the integration constant c 0 = 1 − B(0) is fixed by the initial conditions Φ(0) = 1. Equation (45) then becomes which is the first important result of this study. Together with Eq. (36) we can identify firstly, the relation to earlier obtained solutions for special values of b and k, and secondly derive the general solution for the epidemic time evolution for general values of b and k.

Inverse solution for the general case
In the general case b = k the transcendental Eq. (46) is solved in terms of the realvalued Lambert functions [11]. We discuss below which of the two existing real-valued Lambert functions W 0 (principal) and W −1 (non-principal), respectively, applies in different parameter ranges. For the moment the notation W without an index represents both alternatives.
As Φ ≤ 1 we set Φ = e −x with non-negative values of x to find for Eq. (46) the so-called Lambert equation with The transcendental Eq. (51) has the solution [11] so that where we used the identity e −aW(z) = [W(z)/z] a with a = 1 here. Using r from Eq. (52) we finally find the exact relationship between Φ and ψ as with the positive expression We emphasize that E is nonnegative irrespective of the sign of α. At time τ = 0, E evaluates to E(ψ 0 ) = e −1/α . We had already proven that ψ asymptotically reaches sign(α)∞. This implies E ∞ = lim τ→∞ η[e ψ(1−k)/α) ]. Because α < k and k > 0 in general, the (1 − k)/α has the sign of −α, and E ∞ = 0 for any α = 0. Further, because the derivative of E(ψ) with respect to ψ vanishes only at ψ = ln(−k/b), and thus nowhere for positive k and b, the E(ψ) has no extremum with respect to ψ. For the same reason E has a local maximum with respect to time τ, when ψ exhibits a corresponding minimum in time, which is the case for α ∈ (0, 1). For all other α the E monotonically decreases with time from its initial value E(ψ 0 ). Inserting the solution (55) provides for Eq. (36), dψ/dτ = α − Φ, the nonlinear differential equation with the function E from (56). Making use of the initial condition (33) this readily integrates to We thus arrived at an exact analytical solution of the SIRV equations. We call this an inverse solution, as we have expressed τ in terms of ψ and not vice versa. This is a remarkably compact result, especially, as we know already that ψ(τ), for the relevant case of α ∈ [0, 1], is not a monotonous function, while we should come up with a unique ψ for each τ. The solution to this at first sight apparent contradiction is provided by the two realvalued Lambert functions W 0 (x) and W −1 (x), that have different values over a range of negative x ∈ [−e −1 , 0] values. For positive arguments x, only W 0 (x) is real-valued. The range of application of the two Lambert functions will be discussed in Appendix A. Here we state the main outcome of these considerations.
As long as ψ exhibits a minimum, which is the case for all α ∈ (0, 1), Eq. (58) must be interpreted as or alternatively, in a more symmetric fashion, as where the crossover is located at Φ = α, i.e., where ψ = ψ m has reached its minimum, and the time where this minimum occurs, is given by In the absence of a minimum of ψ, i.e., for α / ∈ (0, 1), there is no crossover, the ψ is monotonous, and one can use Eq. (58) throughout. For α < 0 the argument of the Lambert function is positive, so that one has to use Eq. (58) with W = W 0 .
For α > 1, the minimum of ψ coincides with ψ 0 , hence τ m = 0, so that only the 2nd term in the 2nd case of Eq. (60) survives, again involving the principal W 0 only. The non-principal Lambert function W −1 thus only plays a role in the case α = k − b ∈ (0, 1).

3.5.
Determination of the minimum value ψ m for α ∈ (0, 1) To evaluate τ given by Eq. (60) for the case of α ∈ (0, 1) we need to specify ψ m . We have noted before (see Eq. (36)) that the function ψ attains its minimum value ψ m at Φ m = α, and that a minimum exists for all α ∈ (0, 1). According to the solution (55) we Both, the principal and non-principal, Lambert functions have the same value −1 at the argument (−e −1 ) so that E m = α/e. This yields with Eq. (54) the equation Because ψ m is the value at the minimum, ψ m ≤ ψ 0 automatically holds. This nonlinear equation for ψ m cannot be solved analytically in general. For some special cases such as b = k/2 one can still write down an analytical solution. An approximant for ψ m will be derived below.

Approximated reduction of the exact solution 4.1. Approximate inverse solution τ(ψ)
Here we derive an approximant for the exact inverse solution τ(ψ), that can be later inverted exactly to obtain ψ(τ) in the next subsection. Upon introducing z = −E(x)/α with E(x) given by Eq. (56) the previous inverse solutions (58) and (60) for τ, and also τ m defined by Eq. (61), are of the form where ψ 1 , ψ 2 and µ = 0, or µ = −1 are treated as arbitrary coefficients for the time being, as W µ stands for any of the two Lambert functions, so that τ is of the form (63) for any α / ∈ {0, 1}. Evaluating the required derivative of z with respect to x gives The 1st term can be approximated by unity when x 1 and thus e −x 1, and k/α not too close to unity, i.e., b not too small. The precise range of validity of this approximation will be worked out in Appendix B and section 4.4, where we specify a critical b c below which the current approximation need not be used. For α > 1 and α < 0 one has e −x < e −ψ 0 so that this approximation is applicable for any η 1 when α / ∈ [0, 1]. For α ∈ (0, 1) one has e −x < e −τ m so that the approximation is excellent as long as τ m 1. Under such circumstances, i.e., for b > b c , Eq. (64) is well approximated as Then with the substitution z = we w , corresponding to w = W µ (z), and dz/dw = (1 + w)e w we can calculate the integral (66) in closed form as The minus sign is kept in nominator and denominator as the W µ is typically negative. For α ∈ (0, 1) the ψ initially decays with time, until it reaches its minimum ψ m . At the minimum W µ (z(ψ m )) = −1. Because W µ (z) = 0 only for z = 0, and because E(x) is positive for all η > 0, Eq. (67) has removed the problem with the pole that occurs at z = e −1 in the starting Eq. (63).
It is important to realize that the value for ψ m has completely disappeared within the approximate case. Still, we can use a very similar approximation done here to obtain an explicit expression for ψ m , that might be helpful for the evaluation of the exact inverse solution and the exact ψ m as long as α ∈ (0, 1). If b is not too small, we can approximate the 1 + e ψ m in the Eq. (62) determining ψ m by e ψ m and solve for ψ m analytically. This yields (Fig. 2) Note that this value has the feature W −1 [−E(ψ m )/α] = −1 for η = 0. In practise, for η 1, the ln(1 − η) term can be safely neglected and we shall use while it possible to add this ln(1 − η) correction throughout the rest of this document.

Approximate direct solution ψ(τ)
Next we invert the approximate inverse solution to come up with an approximate direct solution of the SIRV equations. With Eqs. (69)-(71) we have provided approximate expressions for τ m , and τ as function of ψ, for any negative or positive α. Fortunately, we can proceed and invert the relationships without any further assumption to come up with the corresponding, much more convenient, explicit solution ψ(τ). Because we have expressed all SIRV functions in terms of ψ, we then have access to all these quantities in terms of reduced time τ.
We begin with an illustrative example. In Eq. (70) we have provided τ m in terms of b and α. We can invert the relationship to obtain α from τ m and b as follows To prove the equivalence between Eq. (74) and Eq. (70), one has to just insert Eq. (74) into Eq. (70) and to use the fact that the Lambert function W(x) is the inverse function of x(W) = We W . More specifically, consider y(ζ) = ln[−W(ζ)]. Its inverse is given by To derive Eq. (74) we had thus identified ζ = −e −1/α /α and y(ζ) = bτ m . The same procedure can be applied to all expressions from the previous section using different choices for y and ζ.
One more required ingredient is however an expression for ψ in terms of E. It can be readily deduced using the existing assumption that lead to Eq. (73). To derive this Eq. (73) we assumed that 1 + e ψ m ≈ e ψ m . Since the minimum ψ m ≤ ψ at all times, the assumption 1 + e ψ ≈ e ψ is even more appropriate within all remaining times. With this replacement Eq. (56) becomes and this can be solved for ψ to arrive at with E (τ) = E(ψ(τ)) and where we used Eq. (72). We here introduced the symbol E instead of E only to highlight the different argument and to avoid potential confusion. Note that this generalization of Eq. (73) is compatible with the special case ψ m = ψ(τ m ) because E (τ m ) = E(ψ m ) = α/e. As before it is even more convenient to drop a term of O(η) so that ψ(0) coincides with the exact value. Since E (0) = E(ψ 0 ) = e −1/α , the correction results in the simpler Replacing y and ζ, we thus find for α ∈ (0, 1) Note that there is no need to consider two regimes before and after the peak anymore! While the two Lambert functions W 0 and W 1 are very different, they share a common inverse, and Eq. (78) is valid over the whole τ ≥ 0 range. To summarize, upon inserting Eq. (78) into Eq. (77), we end up with a final expression for ψ(τ) for all α ∈ (0, 1), . Upon inserting this term into Eq. (77), one arrives at the explicit approximate solution of the SIRV equations for α / ∈ (0, 1) Not only is this result exactly of the form we obtained for α ∈ (0, 1), it is moreover valid also for the special values of α = 0 and α = 1, as all -at first glance problematic -divergencies have dropped out. From Eq. (80) we can see that ψ ∞ = lim τ→∞ ψ(τ) = sign(α)∞ except for α = 0 or b = k, where ψ ∞ = ψ 0 − 1/k approaches a finite value. The full expression valid for all α is thus in full accord with the exact SIRV solution.

Time-dependency of all remaining SIRV quantities
With ψ(τ) given by Eq. (80) and the corresponding for all values of α at hand it is now straightforward to write down all remaining SIRV quantities, as we have expressed them in terms of ψ above. It is perhaps interesting to note that the approximant shares the most relevant features with the exact solution: ψ(τ m ) = ψ m , ψ(0) = ψ 0 , Φ(0) = 1, and ψ (0) = α − 1 (as required by Eq. (37)) for all α. In the limit of infinitely long times, Φ ∞ = 0, according to Eq. (82) unless b = 0. All remaining SIRV quantities are obtained using ψ(τ) and Φ(τ) from Eqs. (80) and (82) via For α = 0, where we can use ψ(τ) from Eq. (80) with b = k, these expressions solve the SIRV equations (12)- (14) exactly, as can be verified by direct insertion into Eqs. (12)-(15). An alternative proof is provided in Appendix G.1. Otherwise they solve the SIRV equations to within O(η). The version J(τ) given by Eq. (86) ensures that j = dJ/dτ holds exactly. It is a rather tedious exercise to insert the ψ and Φ into Eqs. (84)-(88). In evaluating the limiting values for τ → ∞ one has to carefully consider the qualitatively different regimes α < 0, α = 0, α ∈ (0, 1) and α ≥ 1, as well as k = 0 when α = 0. This can be done but we refrain from writing down all equations for the approximate explicit solution of the SIRV equations. Instead, we provide in Appendix G exact solutions for special cases, and compare the approximate explicit solution with the exact numerical solution for several cases. To provide an example, inserting ψ (80) and Φ (82) into the expressions for S (83) and I (84) yields Next, we focus on the most relevant. measurable features of the SIRV model, that derive from the analytic approximant.

Critical reduced vaccination rate b
As we will demonstrate below, the approximants (80) and (82) capture the exact solution very well (or also exactly for some special cases like α = 0) except for the regime where b stays below a critical b c . The limiting case of b = 0 is known as SIR model and had been treated elsewhere so that the failure of our approximant does not seem to pose a problem. It is however possible to quantify the range of validity of Eqs. (80) and (82). This is done in Appendix B and leads to a critical value b c in terms of η and k b c = 32πkη 2 for the reduced vaccination rate, beyond which the approximant is very accurate. See Fig.  4 for a plot of b c versus k and η. The critical b c decreases with increasing k and decreasing η, while it is quite sensitive to k for values close to unity. At k = 1, Eq. (91) evaluates to (32πη 2 ) 3/5 ≈ 15.9 × η 6/5 and b c is thus roughly proportional to η in that case. For the remaining range of b values below b c we can make use of the known solution of the SIR model, that corresponds to the SIRV model with b = 0. To be specific, the characteristics of time-evolution such as peak time and height of the differential rate, or the final fraction of infected population are well captured by a simple linear interpolation between SIRV values at b = b c and the SIR values. This will be worked out next to conclude with approximants that are valid for any b.

Peak times and peak amplitudes
While S(τ) decreases monotonically with time, both I(τ) and j(τ) = S(τ)I(τ) exhibit a maximum, whose position and amplitude we can calculate from Eqs. (89) and (90). While the general case of arbitrary η is treated in Appendices E amd F, we here limit the analysis to the relevant case of a small initially (and simultaneously) infected fraction η 1 of the population. In this limit we obtain to leading order for all α implying for the differential fraction of newly infected persons The peak time of the maximum in I(τ) is thus so that a peak in I thus exists only if k < 1. For the fraction of currently infected persons at peak time we have Likewise, the peak time of the maximum in j(τ) is at (thin colored lines in Fig. 5) where j achieves the value or equivalently, as the following expression is much more conveniently evaluated at small b (thin colored lines in Fig. 6) A peak in j thus occurs only for k + b < 1, in agreement with our earlier consideration (see Eq. (10)), as we assume small η 1 here. While the peak time τ The maximum rate (98) is not applicable for b < b c below a critical, very small b c , given by Eq. (91). For small b < b c we can make use of the known [3] exact result j SIR max of the SIR model, reproduced in Eq. (A62), and linearly interpolate as (thick colored lines in Fig. 6) One should keep in mind that the maximum of the j(τ) is not located at τ = τ m , which marks the time of the minimum in ψ(τ). And there is also a remaining apparent contradiction. While the initial slope of j(τ) is positive for k + b < 1 − 2η and j(τ) thus going through a peak at a future time in that case, and taking identical values at different times, the ψ(τ) exhibits a minimum for all α ∈ (0, 1). Under such latter conditions, there are at least two times that exhibit the same ψ value. The apparent contradiction finds its explanation in the fact that j(τ) cannot be expressed in terms of ψ(τ) alone, but also involves the derivative of ψ with respect to τ, which is contained in Φ. In other words, both Lambert functions are sometimes (when k − b ∈ (0, 1) and k + b > 1 − 2η) required to describe a j(τ) that is monotonically decreasing, and a single Lambert function is sometimes (when k − b / ∈ [0, 1] and k + b < 1 − 2η) sufficient to describe a j(τ) that exhibits a maximum. Another way to understand this feature is the fact, that the expressions for the SIRV quantities have the same form irrespective the value for α, i.e., irrespective the occurrence of a minimum in ψ(τ).

Total fraction of infected persons
The differential j(τ) given by Eq. (94) can be integrated to obtain the cumulative fraction J(τ) as shown in Appendix D. For J ∞ we thus obtain (thin colored lines in Fig. 7) in terms of the lower incomplete gamma function γ [12]. For the special case of α = 0 (k = b), the Eq. of b such as b = b c , that enters the following Eq. (102). Using exactly the same interpolation approach as before for j max , the J ∞ within the regime of small b < b c is approximated with the help of Eq. (101) by (thick colored lines in Fig. 7) where J SIR ∞ (k) is the known analytical expression for the SIR model, reproduced in Eq. (A61), and b c is given by Eq. (91).

Differential rate
We next provide an analytical approximant for the time-dependent differential rate j(τ), valid for any b, and small η 1. For b exceeding the critical b c , we can just use the expression (94). The comparison with the analytic result is excellent (Fig. 8a-c). For b < b c , on the other hand, we have shown already that the peak time and peak amplitude are well approximated analytically by a linear superposition between the SIRV approximant evaluated at the critical b c , and the analytical SIR expression. The analytical expression for the full time-dependency of j(τ) in this subcritical regime of b < b c is therefore not just a linear superposition of the j(τ) for SIRV and SIR model, because such superposition would not recover the already determined peak time and height. Instead, as we demonstrate with Fig. 8(d-f), and inspired by the earlier observation that the Gauss model [13,14] captures the differential rate very well, the j(τ) is well described for b < b c by the Gaussian with a width w that is determined by [13] where j max and J ∞ are given by Eqs. (100) and (102)

Time scales
Summarizing our analysis here we emphasize that SIRV-pandemic waves in the case 0 < b α < 1 exhibit a clear asymmetry with respect to peak time in reduced and real time. The fraction of infected (I) and recovered (R) persons as well as the daily rate of new infections (j) vary rapidly on the order of the recovery reduced time scale τ r α −1 k −1 , corresponding to recovery real time scale T r 1/µ(0). Alternatively, the fraction of susceptible (S) and vaccinated (V) persons as well as the sum Φ = S + I vary slowly on the order of the much greater vaccination reduced time scale τ v b −1 , corresponding to the vaccination real time scale T v 1/v(0). Thirdly, the cumulative fraction of newly infected persons J(τ) exhibits an asymmetric time structure determined both by τ r and τ v . These behaviors are clearly visible in Fig. 9(b-d) or alternatively Fig. 9(f-h), where we show the time distribution of all SIRV quantities in one plot for different values of b at k = 0.9. The behavior is qualitatively similar but quantitatively different for other values of k. Comparing the asymmetric time distribution of J(τ) with empirical data then should allow the determination of the two parameters b and k.
For comparison we show in Fig. 9(a,d) the SIR-time distributions. Here definitely no enhanced asymmetry occurs. Apart from the absent V(τ) all SIR quantities vary on the same recovery reduced time scale k −1 . Moreover, the SIR S(τ) saturates at the finite value given by S SIR ∞ = 1 − J SIR ∞ , whereas the SIRV S(τ) approaches zero after infinite time.

Comparison of approximate with exact solutions
We have determined all model parameters such as η, k, and b, using currently available public data [15,16] for the population amount, vaccination rate, daily number of newly deceased persons. The parameters as well as the SIRV prediction are collected in Tab. 1 for various countries.
The exact numerical solution of the SIRV model is compared with the approximant for various typical choices of k and b in Fig. 8 country  Table 1: Analysis using data from 18 Mar 2021. For η, k and a we use the current values for the 2nd wave, that started at t II 0 , all from an online resource, Ref.
[16] The begin t v of the vaccination program and the mean daily fraction v of vaccinated population since then we retrieve from Ref. [15], assuming that each person has to be vaccinated twice, and that the vaccination is effective two weeks after the 2nd shot. The remaining quantities are derived from Eqs. (8)  ∞ assuming no vaccinations, (iv) J b=b ∞ assuming ongoing vaccination at the present rate, (v) J b=2b ∞ assuming the vaccination rate had been twice as large. The t 99% denotes the date at which 99% of the final J ∞ has been reached, andJ max = j max a × 10 5 /N is the number of newly infected persons per 100,000 inhabitants within a single day, at peak time. The difference between J b=b ∞ and J b=0 ∞ is the population fraction that profits from the current vaccination program. For all countries η 1, k ∈ [0.7, 1], α ∈ [0.7, 1], b 1 hold. A daily updated, and extended table containing more numbers such as η, v is part of the supplementary material. These three values capture the qualitative behavior for all b. For b b c , vaccination is basically ineffective in reducing the number of infections. For b b c , the vaccination program is highly effective. The crossover is at b = b c , where the reduction becomes significant, and where it depends roughly linear on b. The critical b c as function of k and η is shown in Fig. 4. While b c basically coincides with η for k → 1, at smaller k the critical b c behaves nonlinear with k and η, as shown.
The exact result (solid black lines) for J ∞ versus b for various k is compared with the approximant (colored lines) in Fig. 7. The gap in the colored curves marks the crossover regime, b = b c . For any b, the exact solution is captured by the SIRV approximant. The same is true for remaining quantities such as peak height j max and peak time τ j max , as demonstrated by Figs. 6 and 5.

Application to real data
To make predictions and draw conclusions from the SIRV model about the ongoing pandemic and the vaccination efforts, we have collected current values of infection, vaccination, and recovery rates, as well as population sizes, by making use of existing online resources. We then applied the SIRV model to calculate the time evolution of all the SIRV quantities, including final number of infected persons, or maximum daily number of newly infected persons. Examples are shown by Tab. 1, while a corresponding, daily updated table that includes even more characteristics of the 2nd pandemic wave is part of our supplementary material.
Besides the input parameters of the SIRV model such as the rates the table offers the dimensionless values for k, b, and α = k − b, the criterion ∆ = 1 − 2η − k − b, the date t v marking the beginning of the vaccination program, and various values for the final population fraction that is getting infected before all population has been either recovered or vaccinated. As long as ∆ has positive values, new pandemic waves can occur. Table  1 indicates that apart from Israel this can happen in all countries considered. Only Israel has applied an high enough vaccination rate so that no further Covid-19 waves can occur. The table lists the population fraction J I ∞ that had been infected up to the end of the first pandemic wave, the cumulative fraction J ∞ (t v ) at the time the first person got vaccinated, the hypothetical cumulative fraction J b=0 ∞ assuming there was no vaccination program, the J ∞ = J b=b ∞ using the current value for the mean number of vaccinated persons that gave rise to b, and the hypothetical J b=2b ∞ assuming vaccinations could have been performed at twice the actual speed. In addition, the table shows the maximum number of newly infected persons per day and per 100,000 inhabitants (J max ), and the date t 99% for which 99% of the ultimately infected population fraction has been reached. Note that J b=2b ∞ is not always much smaller than J ∞ because the vaccination program eventually started after the peak time, and because the first wave cumulative fraction J I ∞ sets already a lower limit. Adding the data for all countries of Tab. 1 to Fig. 4, we end up with with Fig. 10, where we have used another colormap for the purpose of this figure. The brightness represents log 10 (b c ), shown as function of k and η. Circles for countries have been placed at positions k and η according to Tab. 1, and the brightness of a filled circle corresponds to the reduced vaccination rate b, again from Tab. 1. If the brightness of a circle exceeds the one of the background, the vaccination rate resides above the critical b c for that country. This is the case for Israel (ISR) and Great Britain (GBR), and Finland (FIN), for example. Circles darker than background, such as for Belgium (BEL), Mexico (MEX) and Argentina (ARG), highlight cases where b < b c .

Summary and conclusions
With the now available vaccination against COVID-19 it is quantitatively explored how vaccination campaigns influence the mathematical modeling of epidemics. For this purpose the well-known susceptible-infectious-recovered/removed (SIR) epidemic model is extended to the fourth compartment V of vaccinated persons and the vaccination rate v(t) that regulates the relation between susceptible and vaccinated persons. The vaccination rate v(t) competes with the infection (a(t)) and recovery (µ(t)) rates in determining the time evolution of epidemics. In order for a pandemic outburst with rising rates of new infections it is required that k + b < 1 − 2η, where k = µ 0 /a 0 and b = v 0 /a 0 denote the initial ratios of the three rates, respectively, and η 1 is the initial fraction of infected persons.
Apparently for the first time we derive analytical solutions for the time-dependence of all relevant quantities Q ∈ [S, I, R, V, j, J] of the SIRV-model where j =J and J are the daily and cumulative fraction of new infections, respectively. As in our earlier analysis of the SIR-model we eliminate one of the time-dependent rates by using the new reduced time-variable τ defined with the infection rate by dτ/dt = a(t). Moreover, we adopt the semi-time case with constant ratios k = µ(t)/a(t) and b = v(t)/a(t) between the infection, recovery and vaccination rates. This assumption still allows us to account for any given time-dependence of the infection rate with the caveat that the recovery and vaccination rate have exactly the same time dependence as the infection rate apart from their different initial values. Exact analytical inverse solutions t(Q) for all relevant quantities Q of the resulting SIRV-model in terms of Lambert functions are derived. The values of the three parameters k, b and η completely determine the reduced time evolution the SIRV-quantities Q(τ).
These inverse solutions can be approximated with high accuracy yielding the explicit reduced time-dependences Q(τ) by inverting the Lambert functions. For a given timedependence of the infection rate a(t) the real time-dependence Q(t) can be inferred. The inversion of the Lambert functions operates well for all ratios b exceeding a small but finite critical value b c . In the range of b ∈ [0, b c ] we interpolate the solution using the exact SIR-solution derived before at b = 0 and the inverted SIRV-solution at b c . This approach is remarkably accurate as the comparison with the exact numerical solutions of the SIRV-equations indicates. The analytical solutions show that SIRV-pandemic waves in the relevant case 0 < b α < 1 exhibit a clear asymmetric distribution in reduced and real time. The fractions of infected (I) and recovered (R) persons as well as the daily rate of new infections (j) vary rapidly on the order of the recovery reduced time scale τ r α −1 k −1 . Alternatively, the fractions of susceptible (S) and vaccinated (V) persons as well as the sum Φ = S + I vary slowly on the order of the much greater vaccination reduced time scale τ v b −1 . This asymmetric SIRV-time behavior is significantly different from the SIR-time behavior, where no time asymmetry occurs, and, apart from the absent V(τ), all SIR quantities vary on the same recovery reduced time scale k −1 .
Clearly our analytical solutions are superior to all numerical ones in the literature as they allow us to identify the main determining parameters of the epidemic waves and to understand the correlations between various monitored observables. And indeed as demonstrated the three parameters b, k and η fully control the evolution of the pandemic wave in reduced time. Also valuable is the use of our exact analytical solutions as benchmark for solutions obtained by solving the SIRV-equations numerically.
The influence of vaccinations on the total cumulative number and the maximum rate of new infections in different countries is calculated by comparing our results with monitored real time Covid-19 data. The reduction in the final cumulative fraction of infected persons and in the maximum daily rate of new infections is quantitatively determined by using the actual pandemic parameters a(0), k and b in different countries. The corresponding numbers for an hypothetical adopted doubled (as compared to the actual one) vaccination rate are also given which allows to quantitively assess the total and maximum casualities caused by the delayed and low-level vaccination coverage in many countries. Moreover, a new criterion is developed that decides on the occurrence of future Covid-19 waves in these countries. Apart from Israel this can happen in all countries considered. Only Israel has applied an high enough vaccination rate so that no further Covid-19 waves can occur.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Range of application for the two Lambert functions
As discussed before in Appendix G of Ref. [11] there is a single real-valued solution W 0 (y) of Lambert's equation for arguments y ≥ 0, referred to as the principal. There are two real valued solutions for y ∈ [−e −1 , 0]: the principal one W 0 ∈ [−1, 0] and the non-principal solution W −1 ≤ −1. For arguments below y < −e −1 only complex-valued solutions exist which are of no interest here because the function Φ is real-valued.
As the function Φ = I + S ∈ [0, 1] we require for the general solution (55) that In Fig. A1 we plot the two functions −αW 0 (−E/α) and −αW −1 (−E/α) entering the constraint (A1). The black line represents the upper limit (A2) and the red line represents the initial E(ψ 0 ) = e −1/α . For α / ∈ [0, 1] the initial value provides a more restrictive upper limit for E, as E is monotonically decreasing in those cases, as we have proven already.
The lower bound of the constraint (A1) is automatically fulfilled for all negative arguments of the Lambert function. The principal Lambert function applies for −1 ≤ W 0 (−E/α) ≤ 0, corresponding to the range which automatically fulfils the upper bound of the constraint (A1) when α < 1 is smaller than unity. The coloured area in the Fig. A1b represents the constraint (53). When α > 1 or α < 0, the constraint (A1) is still fulfilled as long as E ≤ e −1/α , as this value for E inserted into Eq. (A1) gives unity. As we had already shown, ψ monotonically increases (decreases) with time τ for α > 1 (α < 0), E(ψ) varies monotonically with ψ, and E ∞ = 0. The E(ψ) therefore monotonically decreases with time for both α > 1 and α < 0, and thus stays below e −1/α at all times, since E(ψ 0 ) = e −1/α . The constraint (A1) is therefore not only fulfilled automatically for α ∈ (0, 1), but for all α. The only exception from this statement are α = 0 and α = 1, as we have not discussed these special cases here. Likewise, the non-principal Lambert function potentially applies for W −1 (−E/α) < −1, corresponding to the range −αW −1 (−E/α) > α. Together with the right-hand side of the constraint (A1) we find that the non-principal solution potentially applies as solution  Table 2: Critical value b c of the reduced vaccination rate b to be used in the interpolants (100) and (102). Mentioned for comparison: The best fitted value b fit c , the analytic expression (91) for b c that is used throughout this work, and the rough estimate b † c according to Eq. (A6). This table is for η = 10 −6 , the analytic expression (91) works equally well for any η 1.
in the range The coloured area in Fig. A1a represents the constraint (A4). It is not obvious at first glance why there are regions in α-E-space that fulfill both inequalities (A3) and (A4). In those regions only one of the two Lambert functions can lead to the true Φ. The transition between the two Lambert functions is where they meet, i.e., when Φ = α. When α ∈ (0, 1), the Φ at early times involves W −1 . This continues until a point in time where Φ = α is reached. From then on, Φ is determined by the principal solution W 0 .
Upon inspecting the fitted b c for various k and η we find that one has b c 1, J ∞ (k, b c ) η, and a constant proportionality between J ∞ and b 1/3 c . In light of Eq. (A15) this translates into the following equation for b c k b c e 1−k+k ln k bc with some coefficient c = 32π which we determine empirically (Fig. 2) upon comparing the exact solution with the approximant. This Eq. (A7) is solved for b c as follows This final expression for b c is compared with the fitted b c and the above estimate b † c in Tab. 2. The agreement between fitted and analytic b c is excellent for all η and k.

Appendix C Proof of Eq. (71)
Here we prove that W 0 (z(ψ 0 )) = −1/α for all α / ∈ [0, 1), where we recall that z(ψ 0 ) = −e −1/α /α. This identity was used in the derivation of Eq. (71). Making use of the inverse Lambert function W −1 0 (w) = we w , we can take the inverse on both sides which completes the proof, but one has to be careful here. The identity holds only for those α for which W 0 (z) resides on its real-valued regime, i.e., Because the argument −α −1 of the inverse Lambert function is ≥ −1 only for α < 0 and α ≥ 1, the identity does not hold for α ∈ [0, 1). If W 0 is replaced by W −1 , the identity holds for the opposite case of α ∈ (0, 1].

Appendix D Proofs of Eqs. (101) and (102)
With the approximate expression (94) that is useful up to order O(η 2 ) we obtain for corresponding cumulative number fraction with the substitution ξ = e −bτ /b where we have inserted J(0) = η and used the lower incomplete gamma function [12] defined by After infinite time we readily obtain from (A10) where in the last step we used the recurrence formula [12] γ(1 + s, x) = sγ(s, x) − x s e −x for s = k/b and x = 1/b. For the special case of k = b this expression (A12) indeed reproduces Eq. (A53) up to order O(η 2 ), i.e.
For the special case of k = 0, Eq. (A12) simplifies to As already discussed, the Eq. (A12) is useful when b does not exceed b c . This avoids that the latter expression (A13) exceeds unity and thus fails when b is too small. With the asymptotic behavior of the gamma function we obtain from Eq. (A12) for small values of b 1 where we used Stirling's formula Γ(s 1) (2π) 1/2 s s− 1 2 e −s in the last step.

Appendix E Cumulative fraction of infected persons J(τ) for arbitrary η
Starting from our approximants for S(τ) and I(τ) given by Eq. (27), with ψ(τ) from Eq. (26), i.e., the differential rate j(τ) = S(τ)I(τ) (29) is written as Here we are interested in an expression for the integrated j(τ) for arbitrary initial conditions η, contained in ψ 0 = ln[(1 − η)/η], while a result valid for small η 1 we have already provided with Eq. (A10). To this end we calculate with Eqs. (30) and (A17) Substituting ξ = e −bτ /b then yields To be able to calculate this integral we introduce yet another, but simpler integral U(c) parameterized by c, Because the derivative of U(c) with respect to c evaluates to we can express the cumulative fraction (A19) in terms of U(c) as We are left to calculate the integral (A20). Using the identity that all vanish for τ = 0. We have kept the exponential weight containing ψ 0 in front of the coefficient J n (τ) in Eq. (A25) as this allowed us to see that we need to take into account the first m terms of the sum, and the first m coefficients J m to come up with J(τ) up to order O(η m ). The integral in the coefficient can be expressed in terms of the lower incomplete gamma function γ so that we finally have obtained J(τ) as infinite sum, whose summands are given by While this appears as a tractable expression, it cannot be directly evaluated numerically for small b such as b = b c , because e n/b poses a problem. To circumvent this problem, we make use of Stirling's formula and γ(χ n , ∞) = Γ(χ n ). where eventually we re-inserted χ n from Eq. (A31). For small values of η 1 we only take the first term in this sum providing with α = k − b k which agrees exactly with Eq. (A15). It is worthwhile noticing that for very small b 1 the higher order terms in the sum must be taken into account, even at small η, as they make sure that J ∞ stays below unity. We need to calculate J ∞ only down to b = b c , while the form (91) for b c ensures that J ∞ < 1.
At infinite time τ = ∞, we find from Eqs. (A45), (A47), and (A48) as well as I ∞ = S ∞ = j ∞ = 0, R ∞ + V ∞ = 1 and J ∞ = R ∞ . The only nontrivial quantity is thus J ∞ , which we obtain from Eq. (A52) inserted into Eq. (A51) as J ∞ (k, k) = 1 + k ln η 1 + e ψ 0 −(1/k) For small values of k, J ∞ 1 + k ln η which for k = 0 correctly provides J ∞ = 1. All other expressions can also be readily evaluated in this limit with the help of lim k→0 [1 − e −kτ ]/k = τ. For α = k = b = 0, we thus have ψ(τ) = ψ 0 − τ and Φ(τ) = 1 and all above expressions simplify considerably. For example, This completes the analysis of the special case of equal values of b = k, or equivalently, α = 0. The most noteworthy result is the significant reduction in the final cumulative number of new infections with increasing values of k = b as compared to the SI-case with k = b = 0, cf. Eq. (A61).
which agrees exactly with the integrable Eq. (24) in Ref. [3] yielding and for the rate of new infections j and the corresponding cumulative number J j(τ) = SI = de −G dτ , J(τ) = 1 − e −G (A60) To summarize, one finds [3] for the SIR model, the special case of the SIRV model with b = 0,

.1 Alternative inverse solution
For b = 0 one obtains α = k, so that Eq. (36) reads Likewise, the general solution (55) reduces to with the positive expression Inserting the solution (A64) then provides for Eq. (A63) which with the initial condition (33) readily integrates to the inverse exact solution Hence practically all previous general results obtained in Sect. III and IV also hold here with α replaced by k. For values of k ∈ (0, 1), Eq. (A67) yields where the time where the minimum ψ m occurs, is given by The minimum value ψ m for the case of k ∈ (0, 1) is determined by the condition E 0,m = E 0 (ψ m ) = k/e providing with Eq. (A65) Following the Appendix G.3 we can further reduce the inverse solution (A67) written as with z 0 (x) = −E 0 (x)/k = − A(1 + e −x ) k , then provides for Eq. (A71) The substitution z 0 = we w then yields which is still exact. With the substitution y = e −w the solution (A75) can be further rewritten. Integrals of this form have been approximated in Kröger and Schlickeiser [11].
Likewise, the general solution (55) reduces to with the positive expression