Next Article in Journal
Spin Current in BF Theory
Next Article in Special Issue
SIR-PID: A Proportional–Integral–Derivative Controller for COVID-19 Outbreak Containment
Previous Article in Journal
Dynamics of a Homogeneous and Isotropic Space in Pure Cubic f(R) Gravity
Previous Article in Special Issue
Comparison of Main Covid-19 Outbreaks and Interpretation Based on Age Differences
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Institut für Theoretische Physik, Lehrstuhl IV: Weltraum-Und Astrophysik, Ruhr-Universität Bochum, D-44780 Bochum, Germany
2
Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, Leibnizstr. 15, D-24118 Kiel, Germany
3
Polymer Physics, Department of Materials, ETH Zurich, CH-8093 Zurich, Switzerland
*
Authors to whom correspondence should be addressed.
Submission received: 28 March 2021 / Revised: 29 April 2021 / Accepted: 12 May 2021 / Published: 25 May 2021
(This article belongs to the Special Issue Physics Methods in Coronavirus Pandemic Analysis)

Abstract

:
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.

1. Introduction

Statistical physics of vaccination has a rich history in physics research [1]. The model proposed here is an extension of an often employed model in epidemics research, the susceptible-infectious-recovered/removed (SIR) model. While this model can be regarded as a toy model, it has been proven to at least qualitatively capture important aspects of an epidemic and is still heavily used to determine parameters characterizing epidemics or to forecast the amount of required clinical equipment. In this study, vaccinations are taken into account and the corresponding analytic expressions are derived so that those can be used without a computer at hand. Of course, numerical models, including neural networks and artificial intelligence, are convenient frameworks that allow us to take into account many more effects and data, but the goal of the present study is not to add to this area of modeling research but to deliver some insight and analytic approximants for the solution of a well-defined problem, located in the fields of growing dynamical systems and nonlinear differential equations, which to be considered is relevant in the context of Covid-19 research.
There are several variants [2,3,4,5,6,7,8,9,10] of the SIR model, including stochastic variants [11,12,13,14,15,16,17,18,19,20,21,22,23,24] and there is a large amount of related work in the context of the SIR model in the presence of vaccination. Many works deal with vaccination strategies [13,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43], transmission among an interconnected group or population [44] and vaccination behavior by coupling the epidemic spreading with human decisions [45], or policies [46,47] using the SIR model with limited resources [48]. One of the classical variants focuses on an optimal pulse vaccination strategy [49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71], others focus on thresholds and bifurcations in the dynamic epidemic model [72,73,74,75], or the analysis of the delayed SIR model [15,16,22,30,51,61,76,77,78,79]. There are several approaches that investigate the optimal control of the model with time-dependent or nonlinear functions [49,50,55,59,60,62,67,68,69,70,75,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93] and the SIR model had also been applied to networks [34,94,95,96]. While it is comparably easy to set up a network or agent-based [97,98,99] model or to solve a nonlinear problem numerically, pure analytic works are relatively scarce [100,101,102,103,104,105,106,107,108,109,110,111].
In December 2020, the effective mRNA-based Covid-19 vaccine by the companies Pfizer-BioNTech and Moderna became available. This has led to intensive vaccination campaigns in many countries over the world at different speeds. As two shots per person are needed for nearly 100 percent protection as of 10 February 2021, Israel with a 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.
The purpose of the present paper is to analytically and quantitatively investigate, 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. The analysis is based on the SIR epidemic model [112] augmented by the appropriate vaccination rates leading to the susceptible-infectious-recovered/removed-vaccinated (SIRV) epidemic model. This model is a dynamical system for the 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 (Figure 1). Since each vaccine is characterized by a certain efficiency, (effectively) vaccinated persons are defined to be no longer susceptible to being infected. In the case of negligible vaccination, assuming a constant ratio between infection and recovery rate, considerable improvements of the analytical modelling of epidemics with this compartment model have been achieved [109,110]. In what follows, these improvements are of a frequent use.
Application of this earlier study [113] to the monitored second waves in many countries has indicated initial ( t = 0 ) 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 second 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 [110] 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 ( 0 ) / μ ( 0 ) is considerably smaller than the ratio k: the values b [ 5 × 10 4 , 10 2 ] k are expected. Consequently, the difference α = k b is positive and only slightly smaller than k. One may refer to α as the effective ratio as compared to the ratio k.
It is the purpose of this paper to quantitatively determine the influence of this small reduction of the ratio k on the time evolution of the pandemic wave. As in our earlier study [113], the present 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 [114,115,116,117] but they were more concerned with optimizing the control of the epidemics by vaccination with limited resources. The main goal of this study is to derive analytical solutions for the dynamical SIRV. In order to keep the analysis as simple and transparent as possible, such complicating issues as age grouping, vital dynamics and/or spatial spread effects, recently investigated in the literature [118,119] with numerical solutions, are ignored here.

2. 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 ) read:
S ˙ = a ( t ) S I v ( t ) S ,
I ˙ = a ( t ) S I μ ( t ) I ,
R ˙ = μ ( t ) I ,
V ˙ = v ( t ) S ,
where the dot denotes the t-derivative. The SIRV equations are supplemented by semi-time initial conditions (defining t = 0 ),
S ( 0 ) = 1 η , I ( 0 ) = η , R ( 0 ) = V ( 0 ) = 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 (1).
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 N R ( 0 ) , same for V. More formally, if X ˜ denotes one of the SIRV quantities that fulfills initial condition X ˜ ( 0 ) , the time-evolution of X ˜ is given by the time-evolution of X via X ˜ = χ X for X { S , I } and X ˜ = X ˜ ( 0 ) + χ X for X { R , V } with χ = I ˜ ( 0 ) + S ˜ ( 0 ) . The initial conditions (6) are specified by η ˜ = I ˜ ( 0 ) / χ , and the parameters to be introduced next have to be adjusted as well, k ˜ = χ k and b ˜ = χ 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 by
J ˙ = a ( t ) S I ,
and the fraction of vaccinated persons, V ( t ) , are two quantities that can be more easily measured, and are usually reported. Using Equations (2) and (5) the total cumulative number fraction J of persons, that have been infected up to the time t, is related to the SIRV quantities via
J ( t ) = t J ˙ ( ξ ) d ξ = 1 S ( t ) V ( t ) ,
or, equivalently, upon making use of the sum constraint (1), by the sum of currently infected and currently recovered, J ( t ) = I ( t ) + R ( t ) .
Following previous investigations [109,110], as well as to make sure that the SIRV model has any predictive power, the rate a ( t ) is allowed to have an arbitrary time-dependency in contrast to the rates μ ( t ) and v ( t ) considered to share their time-dependency with a ( t ) . This assumption leaves us with the two time-independent model parameters,
k = μ ( t ) a ( t ) , b = v ( t ) a ( t ) .
As it is demonstrated below, 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
τ = 0 t d ξ a ( ξ ) .

2.1. 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,
J ¨ = a ˙ S I + a ( I S ˙ + S I ˙ ) = J ˙ a ˙ / a + a ( S I ) μ v = J ˙ a ˙ / a + a [ ( S I ) k b ] ,
where the argument t is omitted for all functions, and Equation (9) is used. In order for a pandemic outburst with initially rising rates of new infections J ¨ ( 0 ) > 0 to occur, the bracket on the right-hand side of Equation (11) has to be positive at the starting time t = 0 , so that with initially constant rate values a ˙ ( 0 ) = 0 it is required that [ ( S ( 0 ) I ( 0 ) ] k b > 0 . Making use of the first two initial conditions (6), the outburst condition reads:
k + b < 1 2 η .
As k and b have positive values, the condition (12) 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 comment, let us note that Equation (12) demonstrates that in the SIR model with b = 0 (no vaccination), a pandemic does not occur if the parameter k = 1 . The SIR model correctly indicates that epidemic waves end in the case k = 1 . Therefore, the recent criticism [120] about the SIR model is inappropriate and misguided.

2.2. Reduced Time

Let us consider the case when 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 rates have exactly the same time dependence as the infection rate apart from their different initial values. In terms of the reduced time scale τ (10), SIRV Equations (2)–(5) reduce to:
d S d τ = S I b S ,
d I d τ = S I k I ,
d R d τ = k I ,
d V d τ = b S .
Equations (13) and (14) readily yield:
I = b + d ln S d τ ,
S = k + d ln I d τ ,
providing for Equations (15) and (16):
d d τ R + k ln S + b k τ = 0 ,
d d τ V b ln I b k τ = 0 .
Both equations integrate immediately to
R ( τ ) = k ln S ( τ ) + k ln ( 1 η ) b k τ ,
V ( τ ) = b ln I ( τ ) b ln η + b k τ ,
where the integration constants have been determined with the initial conditions (6) holding now at τ = 0 . In terms of the reduced time τ the differential new number of infected persons j ( τ ) = d J ( τ ) / d τ = J ˙ ( t ) / a ( t ) is
j ( τ ) = d J ( τ ) d τ = S ( τ ) I ( τ ) ,
while the corresponding cumulative fraction is
J ( τ ) = τ j ( τ ) d τ = 1 S ( τ ) V ( τ ) = I ( τ ) + R ( τ ) .
At τ = 0 , the cumulative J ( 0 ) = I ( 0 ) + R ( 0 ) = η . In reduced time, Equation (11) corresponds to
1 j d j d τ = d ln j d τ = S ( τ ) I ( τ ) k b .
So far R, V, J, and j were expressed in terms of S and I. An additional relationship between S and I is provided by the sum constraint (1), 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 Equations (17), (18), (21) and (22) provides for the sum constraint (1):
d ln I d τ d ln S d τ + b ln I k ln S = 1 + b k + b ln η k ln ( 1 η ) = 1 + ( b k ) [ 1 + ln η ( 1 η ) ] + k ln η b ln ( 1 η ) .
In the present paper, analytical solutions of Equation (26) are derived for general non-zero and different values of b k . An implicit analytic solution is obtained that expresses the reduced τ in terms of a parameter ψ , while all SIRV-functions including j are expressed in this parameter. Then, a highly accurate analytical approximation for all SIRV-functions as a function of τ is given. These new analytical solutions reduce in the appropriate limit to the earlier [109,110,113] solutions for the non-vaccination case b = 0 . The non-recovery case k = 0 and the case of k = b are considered as special cases.

3. Dynamics of the Epidemics

3.1. Summary of Results

In light of the following rather lengthy derivation of the solution of the SIRV Equations (10)–(16) in reduced time with parameters k and b, and subject to initial condition, I ( 0 ) = 1 S ( 0 ) = η and R ( 0 ) = V ( 0 ) , let us start with the final result. An implicit exact solution τ = τ ( ψ ) parameterized by ψ is derived, while all SIRV quantities can be expressed in ψ as well. Provided the reduced vaccination rate b exceeds a critical b c , for which the explicit expression (92) is provided, it is shown that the explicit solution of the SIRV equation can be written as follows. With the help of the following expression for a function ψ ( τ ) to be defined in the next Subsection,
ψ ( τ ) = ln 1 η η + ( k b ) τ 1 e b τ b ,
one obtains:
S ( τ ) = e b τ 1 + e ψ ( τ ) , I ( τ ) = e b τ 1 + e ψ ( τ ) ,
R ( τ ) = J ( τ ) I ( τ ) , V ( τ ) = 1 S ( τ ) J ( τ ) ,
where the differential j and cumulative fractions J of infected persons are given by
j ( τ ) = S ( τ ) I ( τ ) ,
J ( τ ) = η + 0 τ j ( τ ) d τ .
If b is smaller than the critical b c , provided by Equation (92), 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, these expressions are derived, additional features of the solution, such as the final values at t , are provided. Afterwards, special cases such as the case of k = b are treated, and features of the SIRV equations are discussed which give rise to Equation (27). The solution (27)–(31) of the SIRV model holds for any semipositive values of k and b.

3.2. Two Useful Functions

Let us start the mathematical analysis by introducing the function,
ψ ( τ ) = ln S ( τ ) I ( τ ) ,
implying
I ( τ ) = S ( τ ) e ψ ( τ ) ,
with the usually positive (for initial fractions η < 0.5 of infected persons) initial value,
ψ 0 ψ ( τ = 0 ) = ln 1 η η .
The first derivative of ψ (32) is given by
d ψ d τ = I d S d τ S d I d τ I S = α ( I + S ) ,
where Equations (13) and (14) are used, and the abbreviation
α = k b
is introduced. Consequently, the function
Φ ( τ ) = α d ψ d τ = I ( τ ) + S ( τ )
is always positive and has values Φ [ 0 , 1 ] as I and S are positively valued fractions at all times.
Equation (22) is identical to I ( τ ) = η exp [ k τ ( V ( τ ) / b ) ] , and, since k > 0 and with V residing in the finite interval [ 0 , 1 ] , one finds that, after infinite time, I = I ( τ = ) = 0 . Similarly, Equation (21) 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. At this limit, Equation (37) means 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 , ψ reaches a constant value.
The initial slope of d ψ / d τ evaluates to
d ψ d τ τ = 0 = ( 1 α ) ,
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.,
d Φ d τ = d I d τ + d S d τ = ( b S + k I ) ,
where Equations (13) and (14) are used again.

3.3. Mathematical Analysis

With Equation (33), Equations (37) and (39) yield:
Φ = S ( 1 + e ψ ) , d Φ d τ = S ( b + k e ψ ) .
The combination of Equations (40) then provides:
d Φ d τ = Φ b + k e ψ 1 + e ψ ,
or equivalently,
d ln Φ d τ + A ( ψ ) = 0 ,
with the function
A ( ψ ) = b 1 + e ψ + k 1 + e ψ = d B ( ψ ) d ψ ,
which is written as the derivative of B given by:
B ( ψ ) = b ln ( 1 + e ψ ) k ln ( 1 + e ψ ) = k ψ α ln ( 1 + e ψ ) .
With the function (44), from Equation (37), one obtains for Equation (42) multiplied by ( d ψ / d τ ) = α Φ :
( Φ α ) d ln Φ d τ d ψ d τ d B d ψ = Φ d ln Φ d τ α d ln Φ d τ d B d τ = d Φ d τ α d ln Φ d τ d B d τ = d d τ [ Φ α ln Φ B ] = 0 ,
with the first integral,
c 0 = Φ α ln Φ B ( ψ ) = 1 B ( 0 ) .
Here, the integration constant c 0 = 1 B ( 0 ) is fixed by the initial condition Φ ( 0 ) = 1 . Equation (46) then becomes
Φ α ln Φ = 1 + k ( ψ ψ 0 ) α ln [ η ( 1 + e ψ ) ] ,
which is the first important result of this study. Along with Equation (37), this allows us, first, to identify the relation with earlier obtained solutions for special values of b and k, and, second, to derive the general solution for the epidemic time evolution for general values of b and k. As soon as, according to Equations (47) and (37), the time dependences of Φ ( τ ) and ψ ( τ ) are inferred, one obtains:
S ( τ ) = Φ ( τ ) 1 + e ψ = α d ψ d τ 1 + e ψ ,
I ( τ ) = Φ ( τ ) 1 + e ψ = α d ψ d τ 1 + e ψ ,
according to Equations (33) and (40). Since R ( τ ) and V ( τ ) were already expressed in terms of S ( τ ) and I ( τ ) in Equations (21) and (22), one obtains:
R ( τ ) = b k τ k ln Φ ( τ ) + k ln [ ( 1 η ) ( 1 + e ψ ) ] ,
V ( τ ) = b k τ + b ln Φ ( τ ) b ln [ η ( 1 + e ψ ) ] .

3.4. Inverse Solution for the General Case

In the general case of b k , the transcendental Equation (47) is solved in terms of the Lambert functions [109]. Below, it is discussed which of the two existing real-valued Lambert functions, W 0 (principal) or W 1 (non-principal), has to be applied in the corresponding parameter ranges. Until then, the W functions are used without any index representing both alternatives.
As Φ 1 , Φ = e x is set, with non-negative values of x, in order to find the so-called Lambert equation for Equation (47):
e x = α ( x r ) ,
with
r = 1 + k ( ψ ψ 0 ) α ln [ η ( 1 + e ψ ) ] .
The transcendental Equation (52) has the solution [109]:
x = r + W ( e r α ) ,
so that
Φ = e r W ( e r α ) = α W e r α ,
where the identity e u W ( z ) = [ W ( z ) / z ] u with u = 1 is employed. Using r, Equation (53), one finally finds the exact relationship between Φ and ψ , namely:
Φ = α W E ( ψ ) α ,
with
E ( ψ ) = η ( 1 + e ψ ) e [ 1 + k ( ψ ψ 0 ) ] α 0 .
To emphasize is that E is nonnegative irrespective the sign of α . At time τ = 0 , E evaluates to E ( ψ 0 ) = e 1 / α . It has been already proven above that ψ asymptotically reaches sign ( α ) . This implies E = lim τ η [ e ψ ( 1 k ) / α ) ] . Since α < 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 α values, E monotonically decreases with time from its initial value E ( ψ 0 ) .
Inserting the solution (56) provides for Equation (37), d ψ / d τ = α Φ , the nonlinear differential equation,
d ψ d τ = α 1 + W E ( ψ ) α ,
with the function E (57). Making use of the initial condition (34), this readily integrates to:
τ = 1 α ψ 0 ψ d x 1 + W E ( x ) α .
One thus arrives at an exact analytical solution of the SIRV equations. This to be called an inverse solution, as τ is expressed in terms of ψ and not vice versa. This is a remarkably compact result, especially, since it is already known that ψ ( τ ) , for the relevant case of α [ 0 , 1 ] , is not a monotonous function. However, one has to come up with a unique ψ for each τ .
The solution to this at first sight apparent contradiction is provided by the two real-valued 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 values of x, only W 0 ( x ) is real-valued. The range of application of the two Lambert functions will be discussed in Appendix A. Here, the main outcome of these considerations is stated.
As long as ψ exhibits a minimum, which is the case for all α ( 0 , 1 ) , Equation (59) must be interpreted as:
τ = 1 α ψ 0 ψ d x 1 + W 1 E ( x ) α , ψ ψ m , τ m + 1 α ψ m ψ d x 1 + W 0 E ( x ) α , ψ ψ m ,
or, alternatively, in a more symmetric fashion, as:
τ = τ m + 1 α ψ m ψ d x 1 + W 1 ( z ( x ) ) , ψ ψ m , 1 α ψ m ψ d x 1 + W 0 ( z ( x ) ) , ψ ψ m ,
where the crossover is located at Φ = α , that is, where ψ = ψ m has reached its minimum, and the time, where this minimum occurs, is given by
τ m = 1 α ψ 0 ψ m d x 1 + W 1 E ( x ) α .
In the absence of a minimum of ψ , that is, for α ( 0 , 1 ) , there is no crossover, the ψ is monotonous, and one can use Equation (59) throughout. For α < 0 the argument of the Lambert function is positive, so that one has to use Equation (59) with W = W 0 .
For α > 1 , the minimum of ψ coincides with ψ 0 , hence, τ m = 0 , so that only the second term in the second case of Equation (61) 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 )

In order to evaluate τ given by Equation (61) for the case of α ( 0 , 1 ) , one needs to specify ψ m . As it was noted above, see Equation (37), the function ψ attains its minimum, ψ m at Φ m = α , and a minimum exists for all α ( 0 , 1 ) . According to the solution (56), one finds the minimum value ψ m from Φ m = α = α W ( E m / α ) , where E m = E ( ψ m ) , yielding W ( E m / α ) = 1 . Both, the principal and the non-principal Lambert functions have the same value 1 at 1 / e , so that E m = α / e . Along with Equation (55) this yields the equation that determines ψ m :
e k α ( ψ m ψ 0 ) 1 + e ψ m = η α e 1 ( 1 / α ) .
Since ψ m is the value at the minimum, the inequality ψ m ψ 0 automatically holds. This nonlinear Equation (63) 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 is derived below.

4. Approximated Reduction of the Exact Solution

4.1. Approximate Inverse Solution τ ( ψ )

Here, an approximant is derived for the exact inverse solution τ ( ψ ) that later is inverted exactly in order to obtain ψ ( τ ) in the next Subsection. Upon introducing z = E ( x ) / α with E ( x ) given by Equation (57), the previous inverse solutions (59) and (61) for τ , and also τ m defined by Equation (62), are of the form
τ = 1 α ψ 1 ψ 2 d x 1 + W μ ( E ( x ) α ) = 1 α z ( ψ 1 ) z ( ψ 2 ) d z ( d z / d x ) [ 1 + W μ ( z ) ] ,
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 (64) for any α { 0 , 1 } . Evaluating the required derivative of z with respect to x gives:
d z d x = 1 1 + e x k α z .
The first 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 in Section 4.4, where a critical b c is specified, below which the current approximation needs not be used. For α > 1 and α < 0 , e x < e ψ 0 , so that this approximation is applicable for any η 1 when α [ 0 , 1 ] . For α ( 0 , 1 ) , e x < e τ m , so that the approximation is excellent as long as τ m 1 . Under such circumstances, that is, for b > b c , Equation (65) is well approximated as
d z d x 1 k α z = b α z ,
with α = k b . Hence, Equation (64) is well approximated by
τ 1 b z ( ψ 2 ) z ( ψ 1 ) d z z [ 1 + W μ ( z ) ] .
Then, with the substitution z = w e w , corresponding to w = W μ ( z ) , and d z / d w = ( 1 + w ) e w , one can calculate the integral (67) in closed form:
τ = 1 b W μ ( z ( ψ 2 ) ) W μ ( z ( ψ 1 ) ) d w ( 1 + w ) e w w e w ( 1 + w ) = 1 b W μ ( z ( ψ 2 ) ) W μ ( z ( ψ 1 ) ) d w w = 1 b ln W μ ( z ( ψ 1 ) ) W μ ( z ( ψ 2 ) ) .
The minus sign is kept in the nominator and denominator as W μ is typically negative. For α ( 0 , 1 ) , ψ initially decays with time until reaches its minimum ψ m . At the minimum, W μ ( z ( ψ m ) ) = 1 . Since W μ ( z ) = 0 only for z = 0 and E ( x ) is positive for all η > 0 , Equation (68) removes the problem with the pole at z = 1 / e at starting point, in Equation (64).
Having arrived at Equation (68), one can now identify placeholders ψ 1 , ψ 2 and μ in all three equations (59) (with μ = 0 ), (61), and (62), to write down τ for the various cases as well as τ m .
Let us start with the case of α ( 0 , 1 ) and ψ ψ m . Upon comparing the first line in Equations (61) and (64), ψ 1 = ψ m and ψ 2 = ψ along with μ = 1 for the regime ψ ψ m need to be used. Since α is positive and, thus, z is negative, using Equation (68), this provides immediately:
τ = τ m + 1 b ln W 1 ( z ( ψ m ) ) W 1 ( z ( ψ ) ) ( α ( 0 , 1 ) , ψ ψ m ) .
Since z ( ψ m ) = 1 / e , W μ ( 1 / e ) = 1 , the term ln [ W μ ( 1 / e ) ] = ln ( 1 ) = 0 vanishes for both μ . The nominator in Equation (69) can thus be replaced by unity. Exactly the same procedure, applied to the remaining cases, yields the final expression for the approximate inverse solution when α ( 0 , 1 ) :
τ = τ m 1 b ln W 1 ( E ( ψ ) / α ) , 0 τ τ m , τ m 1 b ln W 0 ( E ( ψ ) / α ) , τ τ m ,
with E ( ψ ) from Equation (57). Similarly, the time τ m can be read off from Equation (68) by choosing ψ 1 = ψ 0 , ψ 2 = ψ m , and μ = 1 . Hence,
τ m = 1 b ln W 1 e 1 / α α = ln ( α ) b ,
where E ( ψ 0 ) = e 1 / α and E ( ψ m ) / α = 1 / e and the identity α W 1 ( z ( ψ 0 ) ) = 1 (see proof in Appendix C) are used to simplify the expression. Such time τ m only exists for α ( 0 , 1 ) .
There is the remaining case of α [ 0 , 1 ] . Equation (59) implies using ψ 1 = ψ 0 , ψ 2 = ψ , and the principal Lambert function ( μ = 0 ), as already discussed. As z ( ψ ) and, thus, W 0 ( z ( ψ ) ) have different signs for α < 0 and α > 1 , both the nominator and denominator in Equation (68) are multiplied by α , to get rid of two different signs for positive and negative values of α . At the same time, because of the identity α W 0 ( z ( ψ 0 ) ) = 1 (see proof in Appendix C), which holds for all α ( 0 , 1 ) , the logarithm of this quantity vanishes, and one arrives at the final expression for the approximate inverse solution valid for all τ 0 and all α [ 0 , 1 ] ,
τ = 1 b ln [ α W 0 ( E ( ψ ) / α ) ] .
Notice that the argument of the logarithm is positive for both α < 0 and α > 1 .
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 term 1 + e ψ m in Equation (63) by e ψ m and solve it for ψ m analytically. This yields (see Figure 2):
ψ m ψ 0 1 α + α ln α α ln ( 1 η ) b .
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 so that
ψ m ψ 0 1 α + α ln α b ,
can be used, while it is possible to add this ln ( 1 η ) correction throughout the rest of this paper.

4.2. Approximate Direct Solution ψ ( τ )

Next, the approximate inverse solution is inverted in order to come up with an approximate direct solution of the SIRV equations. With Equations (70)–(72), the approximate expressions for τ m , and τ as function of ψ are provided for any negative or positive α . Fortunately, one can proceed and invert the relationships without any further assumption to come up with the corresponding, much more convenient, explicit solution ψ ( τ ) . As soon as all SIRV functions are expressed in terms of ψ , all these quantities are then given as function of the reduced time τ .
Let us begin with an illustrative example. In Equation (71), τ m is given in terms of b and α . One can invert the relationship to obtain α from τ m and b, namely:
e 1 / α α = exp b τ m e b τ m .
To prove the equivalence between Equations (71) and (75), one has to just insert Equation (75) into Equation (71) and to use the fact that the Lambert function W ( x ) is the inverse function of x ( W ) = W e W . More specifically, consider y ( ζ ) = ln [ W ( ζ ) ] . Its inverse is given by ζ ( y ) = exp ( y e y ) because ζ ( y ) = exp [ ln ( W ( ζ ) ) e ln ( W ( ζ ) ) ] = W ( ζ ) e W ( ζ ) = ζ . To derive Equation (75) one thus has 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 Equation (74). To derive this Equation (74), one assumes 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, Equation (57) reads:
E ( ψ ) η exp k ψ 0 1 b ψ α ,
and this can be solved for ψ to arrive at
ψ ( τ ) ψ 0 1 α ln ( 1 η ) + α ln [ Ɛ ( τ ) ] b = ψ m α b ln e Ɛ ( τ ) α ,
with Ɛ ( τ ) = E ( ψ ( τ ) ) , and where Equation (73) is used. Here, the symbol Ɛ is introduced instead of E only to highlight the different argument and to avoid potential confusion. Note that this generalization of Equation (74) is compatible with the special case ψ m = ψ ( τ m ) because Ɛ ( τ 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 Ɛ ( 0 ) = E ( ψ 0 ) = e 1 / α , the correction becomes simpler:
ψ ( τ ) = ψ 0 1 + α ln [ Ɛ ( τ ) ] b .
Having expressed ψ in terms of τ and ln [ Ɛ ( τ ) ] , what is left to do is to write down expressions for ln [ Ɛ ( τ ) ] . With such expressions at hand, Equation (78) is the explicit solution of the SIRV equations. The ln [ Ɛ ( τ ) ] depends on the α range.
For α ( 0 , 1 ) , in view of Equation (70), one needs to use y ( ζ ) = b ( τ m τ ) and ζ = z ( ψ ) = E ( ψ ) / α to obtain the explicit ζ ( y ) = exp ( y e y ) , or equivalently, ln [ ζ ( y ) ] = y e y . Replacing y and ζ , one thus finds:
ln [ Ɛ ( τ ) ] = ln ( α ) b ( τ τ m ) e b ( τ τ m )
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 Equation (79) is valid over the whole τ 0 range.
To summarize, upon inserting Equation (79) into Equation (78), one ends up with a final expression for ψ ( τ ) for all α ( 0 , 1 ) ,
ψ ( τ ) = ψ m + α ( τ τ m ) α b 1 e b ( τ τ m ) = ψ 0 + α τ 1 e b τ b ,
where τ m and ψ m were used from Equations (71) and (74) to arrive at the second line. Figure 3 compares the approximant of Equation (80) with the exact solution for three different α values. For the remaining cases of α [ 0 , 1 ] , Equation (72) has to be inverted. Repeating the above procedure, one gets: ln [ Ɛ ( τ ) ] = b τ α 1 exp ( b τ ) for α [ 0 , 1 ] . Upon inserting this term into Equation (78), one arrives at the explicit approximate solution of the SIRV equations for α ( 0 , 1 ) :
ψ ( τ ) = ψ 0 + α τ 1 e b τ b .
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 Equation (81), one 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:
ψ = ψ 0 1 k δ α , 0 + ( 1 δ α , 0 ) sign ( α ) ,
in full accord with the exact SIRV solution.

4.3. Time-Dependency of All Remaining SIRV Quantities

With ψ ( τ ) given by Equation (81) and the corresponding
Φ ( τ ) = α d ψ d τ = e b τ
for all values of α , it is now straightforward to write down all remaining SIRV quantities, as they have been expressed 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 Equation (38)) for all α . In the limit of infinitely long times, Φ = 0 , according to Equation (83) unless b = 0 .
All remaining SIRV quantities to be obtained from Equations (81) and (83) using ψ ( τ ) and Φ ( τ ) :
S ( τ ) = Φ ( τ ) 1 + e ψ ( τ ) ,
I ( τ ) = Φ ( τ ) 1 + e ψ ( τ ) ,
j ( τ ) = S ( τ ) I ( τ ) ,
J ( τ ) = η + 0 τ j ( τ ) d τ ,
R ( τ ) = J ( τ ) I ( τ ) ,
V ( τ ) = 1 S ( τ ) J ( τ ) .
For α = 0 , for which one can use ψ ( τ ) from Equation (81) with b = k , these expressions solve the SIRV Equations (13)–(15) exactly, as can be verified by direct insertion into Equations (13)–(16). An alternative proof is provided in Appendix G.1. Otherwise, they solve the SIRV equations to within O ( η ) . The version J ( τ ) given by Equation (87) ensures that j = d J / d τ holds exactly.
It is a rather tedious exercise to insert the ψ and Φ functions into Equations (85)–(89). 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, in Appendix G, exact solutions for special cases are provided and the approximate explicit solution are compared with the exact numerical solution for several cases. To provide an example, inserting ψ (81) and Φ (83) into the expressions for S (84) and I (85) yields:
S ( τ ) = e b τ 2 1 + tanh α b τ + b ψ 0 + e b τ 1 2 b ,
I ( τ ) = e b τ 2 1 tanh α b τ + b ψ 0 + e b τ 1 2 b .
Next, let us focus on the most relevant, measurable features of the SIRV model, that to be derived from the analytic approximant.

4.4. Critical Reduced Vaccination Rate b c

As it is demonstrated below, the approximants (81) and (83) 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 the approximant obtained does not seem to pose a problem. It is, however, possible to quantify the range of validity of Equations (81) and (83). This is done in Appendix B and leads to a critical value b c in terms of η and k,
b c = 32 π k η 2 3 / 5 exp W 0 6 ( 1 k + k ln k ) 5 ( 32 π k η 2 ) 3 / 5 ,
for the reduced vaccination rate, beyond which the approximant is very accurate; see Figure 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 values close to unity. At k = 1 , Equation (92) 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 the b values below b c one 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 is worked out next to conclude with approximants that are valid for any b.

4.5. Peak Times and Peak Amplitudes

While S ( τ ) decreases monotonically with time, both I ( τ ) and j ( τ ) = S ( τ ) I ( τ ) exhibit a maximum, whose position and amplitude one can calculate from Equations (90) and (91). Whereas the general case of arbitrary η is treated in Appendix E and Appendix F, here, let us limit the analysis to the relevant case of a small initially (and simultaneously) infected fraction η 1 of the population. In this limit, one obtains to leading order for all α :
I ( τ ) η e 1 e b τ b k τ + O ( η 2 ) ,
S ( τ ) e b τ + O ( η ) ,
implying for the differential fraction of newly infected persons,
j ( τ ) η e 1 e b τ b ( k + b ) τ + O ( η 2 ) .
The peak time of the maximum in I ( τ ) is thus
τ max I = ln ( k ) b + O ( η ) ,
so that a peak in I exists only if k < 1 . For the fraction of currently infected persons at peak time, one has:
I max = I ( τ max I ) η e 1 / b k e k / b + O ( η 2 ) .
Likewise, the peak time of the maximum in j ( τ ) is at (thin colored lines in Figure 5)
τ max j = ln ( k + b ) b + O ( η ) ,
where j achieves the value
j max ( k , b ) = η e 1 / b k + b e ( k + b ) / b ,
or, equivalently, as the following expression is much more conveniently evaluated at small b (thin colored lines in Figure 6):
ln j max ( k , b ) = ln ( η ) 1 + 1 k b + k + b b ln ( k + b ) .
Thus, a peak in j occurs only for k + b < 1 , in agreement with the consideration above (see Equation (12)), as small η , η 1 , is assumed here. While the peak time τ max j increases with increasing b in the regime of b < b c , it decays with b for b > b c .
The maximum rate (99) is not applicable for b < b c below a critical, very small b c , given by Equation (92). For small b < b c , one can make use of the known [110] exact result j max SIR of the SIR model, reproduced in Equation (A63), and linearly interpolate as (thick colored lines in Figure 6)
j max ( k , b ) = b b c j max ( k , b c ) + b c b b c j max SIR ( k ) ,
where j max ( k , b c ) and j max SIR ( k ) are given by Equations (99) (evaluated at b = b c ) and (A63), respectively. The analytical expressions for time and amplitude of the daily number of infected persons are two of the more important results of this study. They can be immediately used to rate the effect of vaccination rate for all b, using Equation (99) for b b c and Equation (101) for b < b c , where the critical rate b c is given by Equation (92).
One should keep in mind that the maximum of the j ( τ ) is not located at τ = τ m , which marks the time of the minimum in ψ ( τ ) . 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 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 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 α , that is, irrespective the occurrence of a minimum in ψ ( τ ) .

4.6. Total Fraction of Infected Persons

The differential j ( τ ) given by Equation (95) can be integrated to obtain the cumulative fraction J ( τ ) , as shown in Appendix D. For J , one, thus, obtains (thin colored lines in Figure 7):
J ( k , b ) = η k b b k 1 e 1 b γ k b , 1 b ,
in terms of the lower incomplete gamma function γ [121]. For the special case of α = 0 ( k = b ), Equation (102) agrees with the exact result (A54) up to order O ( η 2 ) . As for j max , one has to distinguish two regimes: (i) the regime of b > b c , where this expression (102) is useful, (ii) the regime of small b < b c , for which one needs, on one hand, the known exact result of the SIR model, J SIR , and on the other, the value for J given by Equation (102), evaluated at b = b c . Since b c is so very small, the direct insertion into (102) is numerically impossible. Therefore, a limiting expression is derived in Equation (A16) of Appendix D, valid for small values of b such as b = b c , that enters Equation (103) just below. 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 Equation (102) by (thick colored lines in Figure 7):
J ( k , b ) = b η b c 2 π k b c e 1 k + k ln k b c + b c b b c J SIR ( k ) ,
where J SIR ( k ) is the known analytical expression [110] for the SIR model, reproduced in Equation (A62), and b c is given by Equation (92).

4.7. Differential Rate

Next, an analytical approximant for the time-dependent differential rate j ( τ ) is provided, valid for any b, and small η , η 1 . For b exceeding the critical b c , one can just use the expression (95). The comparison with the analytic result is excellent, see Figure 8a–c. For b < b c , on the other hand, it was already shown above 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 it is demonstrated in Figure 8d–f and inspired by the earlier observation that the Gauss model [99,122] captures the differential rate very well, the j ( τ ) is well described for b < b c by the Gaussian,
j ( τ ) = j max exp ( τ τ max j ) 2 w 2 ,
with the width w, determined by [99]
w = J π j max ,
where j max and J are given by Equations (101) and (103), respectively.

4.8. Time Scales

Summarizing the above-made analysis, let us emphasize that SIRV-pandemic waves in the case 0 < b α < 1 exhibit a clear asymmetry with respect to the 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 ) . Third, 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 Figure 9b–d or, alternatively, in Figure 9f–h, where the time distribution of all SIRV quantities is shown 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 both the parameters b and k.
For comparison, Figure 9a,d shows 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 as a consequence of vaccinations.

5. Comparison of Approximate with Exact Solutions

All model parameters, such as η , k, and b, were determined using currently available public data [123,124] for the population amount, vaccination rate, daily number of newly deceased persons. The parameters as well as the SIRV prediction are collected in Table 1 for various countries. It should be noted here that there is some significant variation of the experimental data depending on the data source used, and that the determined model parameters are quoted here without error bars for the reason that the differences in the data bases have not yet been fully resolved. On the other hand, a statistical error can be estimated from the variation of the model parameters with time. For this reason, the daily updated model parameters are determined and reported as part of our Supplementary Material.
The exact numerical solution of the SIRV model is compared with the approximant for various typical choices of k and b in Figure 8. To highlight the effect of reduced vaccination rate b, six cases for each k are shown: b = 2 b c , b = 3 b c , b = 5 b c , vanishing b (SIR model), b = b c / 5 , and b = b c / 2 . These six 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 linearly on b. The critical b c as function of k and η is shown in Figure 4. While b c basically coincides with η for k 1 , at smaller k, the critical b c behaves nonlinear with k and η , as it is shown.
In Figure 7, the exact result (solid black lines) for J versus b for various k is compared with the approximant (colored lines). 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 τ max j , as it is demonstrated by Figure 5 and Figure 6.

6. Application to Real Data

To make predictions and to draw conclusions from the SIRV model about the ongoing pandemic and the vaccination efforts, current values of infection, vaccination and recovery rates are collected, as well as population sizes, by making use of existing online resources [123,124]. Then, the SIRV model is applied 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 in Table 1, while a corresponding, daily updated table that includes even more characteristics of the second pandemic wave is part of the Supplementary Material.
Besides the input parameters of the SIRV model, such as the rates, Table 1 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 (ISR) this can happen in all countries considered. Only Israel has applied a high enough vaccination rate so that no further Covid-19 waves can occur. Table 1 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 = 2 b assuming vaccinations could have been performed at twice the actual speed. In addition, Table 1 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 = 2 b 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 Table 1 to Figure 4, one ends up with with Figure 10, where another colormap is used for the specific 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 Table 1, and the brightness of a filled circle corresponds to the reduced vaccination rate b, also taken from Table 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, e.g., for Israel (ISR), Great Britain (GBR), and Finland (FIN). Circles darker than background, such as for Belgium (BEL), Mexico (MEX) and Argentina (ARG) highlight the case when b < b c .

7. 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 time t-dependent 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 to occur, 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, and η 1 is the initial fraction of infected persons.
Apparently for the first time, analytical solutions for the time-dependence of all relevant quantities Q [ S , I , R , V , j , J ] of the SIRV model are derived, where j and J are the daily and cumulative fraction of new infections, respectively. Similar to the earlier analyses [109,110] of the SIR model, one of the time-dependent rates is eliminated by using the new reduced time-variable τ defined with the infection rate by d τ / d t = a ( t ) . Moreover, the semi-time case is adopted with assumed-to-be-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 time-dependence 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 ] the solution is interpolated 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 / a 1 / k . 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 1 / b . 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, the analytical solutions obtained here 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. The use of these exact analytical solutions is also valuable as a 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 the results obtained here 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 a hypothetical adopted doubled (as compared to the actual one) vaccination rate are also given, which allows us to quantitatively assess the total and maximum casualties 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 in 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.

Supplementary Materials

A real-time analysis of Covid-19 data along the lines described in this work is available at https://www.complexfluids.ethz.ch/SIRV.

Author Contributions

Both authors contributed equally to this work. Both authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are openly available in [123,124].

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. [109], 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 ] , for the general solution (56), it is required that
α W E α [ 0 , 1 ] ,
with E > 0 always. While for α < 0 the only real-valued Lambert function is W 0 , for the values α > 0 being of a particular interest here, the argument of the Lambert function is always negative. For those α values, real-valued solutions exist only if
E α e 0.368 α .
Figure A1 shows 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 already proven above.
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
α W 0 E α [ 0 , α ] ,
which automatically fulfils the upper bound of the constraint (A1) when α < 1 . The coloured area in Figure A1b represents the constraint (52). When α > 1 or α < 0 , the constraint (A1) is still fulfilled as long as E e 1 / α , as this value for E, inserted into Equation (A1), gives unity. As already shown above, ψ 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 , but these special cases are not discussed here.
Figure A1. Values of the (left) non-principal and (right) principal Lambert functions determining Φ in Equation (56). The solid black line represents the upper limit, E = α / e = 0.368 α (for all α > 0 ), where Φ = α . Only for the coloured region below this line real-valued solutions of Equation (56) exist. The red line represents the initial E ( ψ 0 ) , which serves as a more restrictive upper limit for α > 1 . In their coloured areas, the corresponding Lambert functions apply, while the white regions cannot be reached. Points residing in ( α , E )-space that are shared by both Lambert functions are visited at different times, as exemplarily shown by the white arrows for α = 0.4 . At earlier times the non-principal Lambert function W 1 describes Φ , while W 0 overtakes at later times. The crossover occurs at Φ = α . At this moment, both Lambert function exhibit exactly the same value. For negative α < 0 , the accessible E span the huge range E [ 0 , e 1 / α ] and are therefore not shown. For such α < 0 , the colored region exists only for W 0 .
Figure A1. Values of the (left) non-principal and (right) principal Lambert functions determining Φ in Equation (56). The solid black line represents the upper limit, E = α / e = 0.368 α (for all α > 0 ), where Φ = α . Only for the coloured region below this line real-valued solutions of Equation (56) exist. The red line represents the initial E ( ψ 0 ) , which serves as a more restrictive upper limit for α > 1 . In their coloured areas, the corresponding Lambert functions apply, while the white regions cannot be reached. Points residing in ( α , E )-space that are shared by both Lambert functions are visited at different times, as exemplarily shown by the white arrows for α = 0.4 . At earlier times the non-principal Lambert function W 1 describes Φ , while W 0 overtakes at later times. The crossover occurs at Φ = α . At this moment, both Lambert function exhibit exactly the same value. For negative α < 0 , the accessible E span the huge range E [ 0 , e 1 / α ] and are therefore not shown. For such α < 0 , the colored region exists only for W 0 .
Physics 03 00028 g0a1
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), one finds that the non-principal solution potentially applies as solution in the range
α W 1 E α [ α , 1 ] .
The coloured area in Figure A1a represents the constraint (A4). It is not obvious at first glance why there are regions in the α 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 ) , Φ 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 .

Appendix B. The Critical Vaccination Rate b c

It is straightforward to estimate the value b c of the parameter b, below which the approximations (81) and (83) break down. This occurs when the approximation in Equation (66), based on the assumption of values of ψ 1 , is no longer valid, i.e., when the minimum ψ m of Equation (74) exceeds a certain value of order unity. Consequently, the approximation tends to break down if the following equation is fulfilled for b c ,
1 = ψ 0 1 ( k b c ) + ( k b c ) ln ( k b c ) b c .
This Equation (A5) can be cast into the form of the Lambert Equation (52) with the solution
b c = k e 2 ψ 0 exp W 0 ( 1 + k k ψ 0 ) e 2 ψ 0
in terms of the principal Lambert function W 0 . The solution is denoted here by b c , as the b c to be used is slightly different, as it is shown below. The non-principal Lambert function would produce large b c values for which the approximants do not require any special treatment. Using ψ 0 = ln [ ( 1 η ) / η ] , the expression (A6) may be rewritten further. It is important to note that the approximants for the regime b < b c cannot be evaluated anymore for values well below b c , especially for k very close to unity. In Table A1, therefore, the best values for b c , that were obtained upon fitting the exact solution with the obtained approximant, are collected and compared with those of expression (A6) as well as with another formula to be derived next, that will actually be used as it captures the best b c with much higher accuracy than b c from Equation (A6).
Table A1. Critical value b c of the reduced vaccination rate b to be used in the interpolants (101) and (103). Mentioned for comparison are: the best fitted value b c fit , the analytic expression (92) for b c that is used throughout this paper and the rough estimate b c according to Equation (A6). The value η = 10 6 is used for this table; the analytic expression (92) works equally well for any η 1 .
Table A1. Critical value b c of the reduced vaccination rate b to be used in the interpolants (101) and (103). Mentioned for comparison are: the best fitted value b c fit , the analytic expression (92) for b c that is used throughout this paper and the rough estimate b c according to Equation (A6). The value η = 10 6 is used for this table; the analytic expression (92) works equally well for any η 1 .
k log 10 ( b c fit ) log 10 ( b c ) log 10 ( b c )
0.100 1.18 1.19 1.18
0.300 1.44 1.44 1.53
0.500 1.75 1.74 1.88
0.650 2.06 2.05 2.25
0.800 2.51 2.50 2.77
0.850 2.74 2.73 3.03
0.900 3.05 3.04 3.39
0.950 3.56 3.57 4.00
0.980 4.21 4.23 4.80
0.990 4.69 4.70 5.41
0.995 5.13 5.13 6.01
Upon inspecting the fitted b c , b c fit , values for various k and η , one finds that b c 1 , J ( k , b c ) η , and there is a constant proportionality between J and b c 1 / 3 . In light of Equation (A16), this translates into the following equation for b c :
k b c exp 1 k + k ln k b c = b c 1 / 3 η C ,
with some coefficient C = 32 π , which is determined empirically (Figure 2) upon comparing the exact solution with the approximant. Equation (A7) is solved for b c as follows:
b c = 32 π k η 2 3 / 5 exp W 0 6 ( 1 k + k ln k ) 5 ( 32 π k η 2 ) 3 / 5 .
This final expression for b c is compared with the b c fit values, and the above estimate b c in Table A1. The agreement between fitted and analytic b c is found to be excellent for all η and k (Table 1 shows the comparison only for a selected, representative η value).

Appendix C. Proof of Equation (72)

Here, it is proved that W 0 ( z ( ψ 0 ) ) = 1 / α for all α [ 0 , 1 ) , where to recall is that z ( ψ 0 ) = e 1 / α / α . This identity was used in the derivation of Equation (72). Making use of the inverse Lambert function W 0 1 ( w ) = w e w , one can take the inverse on both sides:
z ( ψ 0 ) = e 1 / α α = W 0 1 1 α = 1 α e 1 / α ,
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., W 0 [ 1 , ] for z [ e 1 , ] . 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 Equations (102) and (103)

With the approximate expression (95), that is useful up to order O ( η 2 ) , one obtains for corresponding cumulative number fraction, with the substitution ξ = e b τ / b :
J ( τ ) = J ( 0 ) + 0 τ j ( τ ) d τ = η 1 + b k b e 1 b e b τ b 1 b ξ k b e ξ d ξ = η 1 + b k b e 1 b γ 1 + k b , 1 b γ 1 + k b , e b τ b ,
where J ( 0 ) = η is inserted and the lower incomplete gamma function [121] defined by
γ ( s , x ) = 0 x ξ s 1 e ξ d ξ
is used. After infinite time, from (A10), one readily obtains:
J ( k , b ) = η 1 + b k b e 1 b γ 1 + k b , 1 b = η k b k b 1 e 1 b γ k b , 1 b ,
where, in the last step, the recurrence formula [121] γ ( 1 + s , x ) = s γ ( s , x ) x s e x for s = k / b and x = 1 / b is used. For the special case of b = k , the obtained expression (A12) indeed reproduces Equation (A54) of Appendix G.1 (where we treat this special case of b = k in detail) up to order O ( η 2 ) , i.e.,
J ( k , k ) = η k ( e 1 / k 1 ) ,
where γ ( 1 , x ) = 1 e x is used. For the special case of k = 0 , Equation (A12) simplifies to
J ( 0 , b ) = e 1 / b η .
As already discussed, Equation (A12) is useful when b does not exceed b c . This avoids that the latter expression (A14) exceeds unity and, thus, fails when b is too small.
With the asymptotic behavior of the gamma function,
γ ( s , x 1 ) Γ ( s ) x s 1 e x Γ ( s ) ,
from Equation (A12), one obtains for small values of b 1 :
J ( k , b 1 ) η k b k b 1 e 1 b Γ k b η 2 π k b exp 1 k + k ln k b ,
where the Stirling’s formula Γ ( s 1 ) ( 2 π ) 1 / 2 s s 1 2 e s is used in the last step.

Appendix E. Cumulative Fraction of Infected Persons J(τ) for Arbitrary η

Starting from the found approximants for S ( τ ) and I ( τ ) given by Equation (28), with ψ ( τ ) from Equation (27), i.e.,
ψ ( τ ) = ψ 0 + α τ 1 e b τ b , α = k b ,
the differential rate j ( τ ) = S ( τ ) I ( τ ) (30) is written as
j ( τ ) = e 2 b τ 4 cosh 2 ψ 2 = e 2 b τ ψ ( τ ) ( 1 + e ψ ) 2 = e ψ ( τ ) 2 b τ ( 1 + e ψ ) 2 .
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 was already provided with Equation (A10). To this end, using Equations (31) and (A18) one gets:
J ( τ ) = η + 0 τ d τ j ( τ ) = η + 0 τ d τ e 2 b τ ψ ( τ ) ( 1 + e ψ ( τ ) ) 2 = η + 0 τ d τ e [ ψ 0 + ( 2 b + α ) τ 1 e b τ b ] ( 1 + e [ ψ 0 + α τ 1 e b τ b ] ) 2 .
Substituting ξ = e b τ / b then yields:
J ( τ ) = η + e b τ b 1 b d ξ ( b ξ ) α b + 1 e 1 b ψ 0 ξ [ 1 + ( b ξ ) α b e 1 b ψ 0 ξ ] 2 .
To be able to calculate this integral, let us introduce yet another but simpler integral U ( c ) , parameterized by c:
U ( c ) = c η + b e b τ b 1 b d ξ 1 + e 1 b ψ 0 ( b ξ ) α b e c ξ .
Because the derivative of U ( c ) with respect to c evaluates to
d U ( c ) d c = η + e b τ b 1 b d ξ e 1 b ψ 0 ( b ξ ) α b + 1 e c ξ ( 1 + e 1 b ψ 0 ( b ξ ) α b e c ξ ) 2 ,
the cumulative fraction (A20) can be expressed in terms of U ( c ) as
J ( τ ) = d U d c c = 1 .
What is left to do is to calculate the integral (A21). Using the identity
1 1 + e x = n = 0 ( 1 ) n e n x ,
applied to Equation (A21), provides:
U ( c ) = c η + n = 0 ( 1 ) n e n b n ψ 0 b n α b + 1 e b τ b 1 b d ξ ξ n α b e c n ξ = c η + 1 e b τ + n = 1 ( 1 ) n e n b n ψ 0 b n n α b + 1 n e b τ b n b d z z n α b e c z ,
where z = n ξ is substituted in the last step. Equation (A25) readily yields, upon replacing τ 0 , and, using the binomials B m , n = m ! / [ n ! ( m n ) ! ] ,
J ( τ ) = d U ( c ) d c c = 1 = η + n = 1 e n ψ 0 J n ( τ ) = η + n = 1 η 1 η n J n ( τ ) = η + m = 1 n = 1 m B m 1 , n 1 J n ( τ ) η m ,
with the coefficients,
J n ( τ ) = ( 1 ) n + 1 e n b b n n α b + 1 n e b τ b n b d z z n α b + 1 e z ,
that all vanish for τ = 0 . The exponential weight containing ψ 0 in front of the coefficient J n ( τ ) in Equation (A26) is kept, as this allowed us to see that one needs 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, finally, J ( τ ) is obtained as an infinite sum, with summands given by
J n ( τ ) = ( 1 ) n + 1 e n b b n n α b + 1 γ n α b + 2 , n b γ n α b + 2 , n e b τ b .
To confirm that the obtained general expression reduces to (A10) already derived for small η 1 , let us evaluate Equation (A26) to first order in η . With B 0 , 0 = 1 , Equation (A26) implies:
J ( τ ) η + J 1 ( τ ) η = η + η e 1 b b α b + 1 γ α b + 2 , 1 b γ α b + 2 , e b τ b = η 1 + e 1 b b k b γ 1 + k b , 1 b γ 1 + k b , e b τ b ,
where J 1 is taken from Equation (A28), and α is replaced by k b . The obtained Equation (A29) is indeed identical to Equation (A10).
One can rewrite J ( τ ) further in a way that helps calculating J = lim τ J ( τ ) using the recurrence formula for the lower incomplete gamma function,
γ ( a + 1 , x ) = a γ ( a , x ) x a e x
Using this recurrence in Equation (A29), one ends up with
J ( τ ) = 1 e b τ + e ψ 0 + k τ 1 e b τ b + n = 1 ( 1 ) n + 1 χ n e n b n ψ 0 b n χ n × γ χ n , n b γ χ n , n e b τ b ,
where the abbreviation,
χ n = 1 + n α b ,
is used only to shorten the expression. Note that χ 1 = k / b . To derive expression (A31), the following identities,
n = 1 ( 1 ) n e n x = 1 1 + e x
as well as
n = 1 ( 1 ) n + 1 e n ψ 0 = n = 1 ( 1 ) n + 1 e n ψ 0 = 1 1 + e ψ 0 = η
are used. After infinite time the first term of Equation (A31) vanishes, and γ ( χ n , 0 ) = 0 can be used. Hence,
J = n = 1 ( 1 ) n + 1 χ n e n b n ψ 0 b n χ n γ χ n , n b .
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, Stirling’s formula and γ ( χ n , ) = Γ ( χ n ) is used. For small values of b, one then obtains:
J b 1 n = 1 ( 1 ) n + 1 χ n e n b n ψ 0 b n χ n Γ ( χ n ) = n = 1 ( 1 ) n + 1 χ n e n b n ψ 0 n α b b n χ n Γ n α b 2 π b α n = 1 ( 1 ) n + 1 χ n α n α b e n [ ψ 0 + α 1 b ] n ,
where, eventually, χ n from Equation (A32) is re-inserted. For small values of η , η 1 , one takes only the first term in this sum providing with α = k b k :
J b 1 , η 1 η 2 π b α k b α k b 1 e 1 α b = η 2 π k b exp 1 k + k ln k b ,
which agrees exactly with Equation (A16). It is worthwhile noticing that for very small b: b 1 the higher order terms in the sum must be taken into account, even at small η , as this ensures that J stays below unity. One needs to calculate J only down to b = b c , while the form (92) for b c ensures that J < 1 .

Appendix F. Peak Time and Amplitude for bbc and Arbitrary η

The approximant j ( τ ) , valid for b b c and given by Equation (A18), can be written as:
j ( τ ) = e A b 2 F ( x ) , F ( x ) = x 2 + p e x ( 1 + e A x p e x ) 2 ,
with p = α / b , e A = e ψ 0 e 1 / b b p = η e 1 / b b p / ( 1 η ) , and x = e b τ / b , so that F ( x ) contains the dependency on time τ via x. Note that x is positive at all times, because b is positive. As soon as x varies monotonically with τ , to find the peak time τ max j and peak height j max , one needs to determine the position and value of F ( x ) at its maximum, provided such maximum exists. For the derivative of F ( x ) with respect to x, one has:
F ( x ) = e x x 1 + p [ ( 2 + p x ) e x + e A x p ( 2 p + x ) [ e 2 + e A x p ] 3 .
Hence, the derivative vanishes at x = x max , where x max solves the highly nonlinear equation,
e x max A x max p = p x max 2 p x max + 2
or equivalently and more suitable for any numerical implementation,
x max A p ln ( x max ) = ln p x max 2 p x max + 2 .
Inserting Equation (A41) into Equation (A38) then provides for the differential rate at peak time:
j max = j ( x max ) = b x max 2 2 1 2 p x max 2 .
This is an expression valid for arbitrary η and b b c in terms of the solution x max of Equation (A41), that can only be obtained numerically. For the special case of small b k , one can proceed analytically that would provide an approximate solution for x max .

Special Case of bc ≤ b ≪ k

For this case of small b k , the above p-variable (see (A38)) becomes large and is well approximated by p k / b . At the same time, the right hand side of Equation (A41) vanishes. The resulting equation for x max ,
x max p e x max = e A ,
is solved with the help of Lambert’s principal function as:
x max = α b W 0 e A / p p = α b W 0 e ( 1 b ψ 0 ) α α ,
where p is replaced and properties of the Lambert function are used. The corresponding j max is still given by Equation (A42). Because ψ 0 is of order unity, and if b ψ 0 1 holds as well, the argument of W 0 simplifies for b k to e 1 / k / k . This latter argument equals 1 / e to within 2.5% for all k [ 0.8 , 1 ] , so that x max α / b = p , and equivalently, τ max ln ( α ) / b can be used under such circumstances.

Appendix G. Exact Solutions for Special Cases

Appendix G.1. The Equal Value Case b = k Corresponding to α = 0

The general analysis above can also be used for a number of special cases to be investigated in this and the next Appendices. Let us start with the special case α = 0 , so that the general Equations (37) and (47) simplify to
Φ = d ψ d τ = 1 + k ( ψ ψ 0 ) .
With the initial condition ψ ( 0 ) = ψ 0 , Equation (A45) immediately integrates to
ψ ( τ ) = ψ 0 + e k τ 1 k ,
implying
Φ ( τ ) = I ( τ ) + S ( τ ) = d ψ d τ = e k τ .
This proves that, for special case of α = 0 , Equation (80) is actually exact and not an approximant. With Equations (A46) and (A47), one readily obtains for Equations (48) and (49):
I ( τ ) = e k τ 1 + e ψ = e k τ 2 1 tanh ψ 2 = e k τ 1 + exp [ ψ 0 1 e k τ k ] ,
S ( τ ) = e k τ 1 + e ψ = e k τ 2 1 + tanh ψ 2 = e k τ 1 + exp [ 1 e k τ k ψ 0 ] .
Consequently, the rate of new infections is then given by
j ( τ ) = S ( τ ) I ( τ ) = e 2 k τ + ψ ( τ ) ( 1 + e ψ ) 2 = 1 4 e 2 k τ cosh 2 1 e k τ 2 k ψ 0 2 .
Using Equation (51) and Φ from Equation (A47), V ( τ ) is written in terms of ψ as:
V ( τ ) = k ln [ η ( 1 + e ψ ) ] .
With the help of the sum constraint (1), the remaining R ( τ ) is given by
R ( τ ) = 1 e k τ + k ln [ η ( 1 + e ψ ) ] .
With J ( τ ) = I ( τ ) + R ( τ ) , there is no need to integrate j ( τ ) in order to come up with a final expression for J ( τ ) using I ( τ ) and R ( τ ) from Equations (A48) and (A52).
At infinite time τ = , from Equations (A46), (A48), and (A49), one finds:
ψ ( k , k ) = ψ 0 1 k ,
as well as I = S = j = 0 , R + V = 1 and J = R . The only nontrivial quantity is, thus, J , which one obtains from Equation (A53), inserted into Equation (A52):
J ( k , k ) = 1 + k ln η 1 + e ψ 0 ( 1 / k ) = k ln η + k ψ 0 + k ln ( 1 + e 1 k ψ 0 ) = k ln ( 1 η ) + k ln ( 1 + e 1 k ψ 0 ) .
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 , one, thus, gets: ψ ( τ ) = ψ 0 τ and Φ ( τ ) = 1 , and all above expressions simplify considerably. For example,
J k = 0 ( τ ) = 1 2 1 + tanh τ ψ 0 2 .
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. Equation (A62).

Appendix G.2. SIR-Case b = 0, k > 0

For no-vaccination campaigns ( b = 0 ), the SIRV model reduces to the SIR model analyzed before [110,113]. For b = 0 , implying α = k 0 , Equation (47) simplifies to:
Φ k ln Φ + k ln ( 1 + e ψ ) = 1 k ln ( 1 η ) ,
or, with Φ = I + S and Equation (48), to:
I + S k ln S = 1 k ln ( 1 η ) .
For b = 0 , Equation (17) simplifies to I = d ln S / d τ , so that Equation (A57) becomes:
d ln S d τ + S k ln S = 1 k ln ( 1 η ) .
In terms of the positively-valued function G = ln ( S ) , Equation (A58) reads:
d G d τ = 1 e G k G k ln ( 1 η ) ,
which agrees exactly with the integrable Equation (24) in Ref. [110] yielding:
I ( τ ) = d G d τ , S ( τ ) = e G ,
and for the rate of new infections j and the corresponding cumulative number J:
j ( τ ) = S I = d e G d τ , J ( τ ) = 1 e G .
To summarize, one finds [110] for the SIR model, the special case of the SIRV model with b = 0 ,
J SIR = 1 + k W 0 1 η k e 1 k ,
j max SIR = k 2 4 1 + W 1 2 ( 1 η ) k e 1 k + 1 2 1 .

Alternative Inverse Solution

For b = 0 , one obtains α = k , so that Equation (37) reads:
d ψ d τ = k Φ .
Likewise, the general solution (56) reduces to
Φ = k W E 0 ( ψ ) k
with the positive expression
E 0 ( ψ ) = η ( 1 + e ψ ) e ψ 0 1 k = ( 1 η ) e 1 k ( 1 + e ψ ) 0 ,
Inserting the solution (A65) then provides for Equation (A64):
d ψ d τ = k 1 + W E 0 ( ψ ) k .
which, with the initial condition (34), readily integrates to the inverse exact solution:
τ = 1 k ψ 0 ψ d x 1 + W E 0 ( x ) k .
Hence, practically all previous general results obtained in Section 3 and Section 4 also hold here with α replaced by k. For values of k ( 0 , 1 ) , Equation (A68) yields:
τ = τ m + 1 k ψ m ψ d x 1 + W 1 ( E 0 ( x ) k ) , ψ ψ m , 1 k ψ m ψ d x 1 + W 0 ( E 0 ( x ) k ) , ψ ψ m ,
where the time where the minimum ψ m occurs, is given by
τ m = 1 k ψ 0 ψ m d x 1 + W 1 E 0 ( x ) k .
The minimum value ψ m for the case of k ( 0 , 1 ) is determined by the condition E 0 , m = E 0 ( ψ m ) = k / e along with Equation (A66):
ψ m = ln k e 1 k 1 1 η 1 .
Following Appendix G.3, one can further reduce the inverse solution (A68) written as:
τ = 1 k ψ 1 ψ 2 d x 1 + W μ ( E 0 ( x ) k ) = 1 k z 0 ( ψ 1 ) z 0 ( ψ 2 ) d z 0 ( d z 0 / d x ) [ 1 + W μ ( z 0 ) ]
with
z 0 ( x ) = E 0 ( x ) / k = A ( 1 + e x ) k , A = ( 1 η ) e 1 k .
Evaluating
d z 0 d x = A e x k = A + k z 0 k
then provides for Equation (A72):
τ = z 0 ( ψ 1 ) z 0 ( ψ 2 ) d z 0 ( A + k z 0 ) [ 1 + W μ ( z 0 ) ] .
The substitution z 0 = w e w then yields:
τ = W μ ( z 0 ( ψ 2 ) ) W μ ( z 0 ( ψ 1 ) ) d w k w + A e w ,
as the exact solution of the SIR model considered in the current Appendix G.2. With the substitution y = e w the solution (A76) can be further rewritten. Integrals of this form have been approximated in Ref. [109].

Appendix G.3. SIV-Case b > 0, k = 0

In the case of a negligible recovery rate k = 0 , implying α = b , Equation (47) simplifies to:
Φ + b ln Φ b ln ( 1 + e ψ ) = Φ + b ln Φ 1 + e ψ = 1 + b ln η ,
or with Φ = I 0 + S 0 (the index 0 indicates that the limit k = 0 is considered here) and Equation (49):
I 0 + S 0 + b ln I 0 = 1 + b ln η b .

Appendix G.3.1. Symmetry Argument

By a simple symmetry argument, the solution of Equation (A78) can be expressed in terms of the SIR-function [109] G obeying Equation (A59). Setting
I 0 = S , S 0 = I , b = k , η b = 1 η ,
Equation (A78) reads:
I + S k ln S = 1 k ln ( 1 η ) ,
which is identical to Equation (A57). Therefore, one can use the SIR-function G ( τ ) , obeying Equation (A59), but now with negative values of k, to obtain for the SIV solutions,
S 0 ( τ ) = d G d τ , I 0 ( τ ) = d G ( τ ) d τ ,
j ( τ ) = d e G d τ , J ( τ ) = 1 e G .
If to follow this approach, one needs to calculate the function G for negative values of k.

Appendix G.3.2. Alternative Inverse Solution

For k = 0 , one obtains α = b , so that Equation (37) reads:
d ψ d τ = ( b + Φ ) .
Likewise, the general solution (56) reduces to
Φ = b W 0 E b ( ψ ) b
with the positive expression
E b ( ψ ) = η e 1 b ( 1 + e ψ ) 0 .
As its argument is positive, it has to be the principal Lambert function in Equation (A83). Inserting the solution (A83) then provides for Equation (A82):
d ψ d τ = b 1 + W 0 E b ( ψ ) b ,
which, along with the initial condition (34), readily integrates to the inverse exact solution:
b τ = ψ 0 ψ d x 1 + W 0 E b ( x ) b .
Here also, nearly all previous general results obtained in Section 3 and Section 4 hold with α replaced by b . With
z b ( x ) = E b ( x ) / b = 1 + e x D , D = b η e 1 b ,
implying
x = ln ( D z b 1 ) , d x d z b = D D z b 1 ,
one finds for Equation (A86):
b τ = D 1 η D 1 + e ψ D d z b ( D z b 1 ) [ 1 + W 0 ( z b ) ] .
The substitution z b = w e w then yields:
b τ = D W 0 ( 1 η D ) W 0 ( 1 + e ψ D ) d w D w e w ,
which is an exact expression for the reduced time τ of the SIV model considered in the current Appendix G.3.

References

  1. Wang, Z.; Bausch, C.T.; Bhattacharyya, S.; d’Onofrio, A.; Manfredi, P.; Perc, M.; Perra, N.; Salathe, M.; Zhao, D.W. Statistical physics of vaccination. Phys. Rep. 2016, 664, 1–113. [Google Scholar] [CrossRef]
  2. Cadoni, M.; Gaeta, G. Size and timescale of epidemics in the SIR framework. Phys. D 2020, 411, 132626. [Google Scholar] [CrossRef] [PubMed]
  3. Chekroun, A.; Kuniya, T. Global threshold dynamics of aninfection age-structured SIR epidemic model with diffusion under the Dirichlet boundary condition. J. Differ. Equ. 2020, 269, 117–148. [Google Scholar] [CrossRef]
  4. Imron, C.; Hariyanto; Yunus, M.; Surjanto, S.D.; Dewi, N.A.C. Stability and persistence analysis on the epidemic model multi-region multi-patches. J. Phys. Conf. Ser. 2019, 1218, 012035. [Google Scholar] [CrossRef]
  5. Karaji, P.T.; Nyamoradi, N. Analysis of a fractional SIR model with general incidence function. Appl. Math. Lett. 2020, 108, 106499. [Google Scholar] [CrossRef]
  6. Mohamadou, Y.; Halidou, A.; Kapen, P.T. A review of mathematical modeling, artificial intelligence and datasets used in the study, prediction and management of COVID-19. Appl. Intell. 2020. [Google Scholar] [CrossRef]
  7. Samanta, S.; Sahoo, B.; Das, B. Dynamics of an epidemic system with prey herd behavior and alternative resource to predator. J. Phys. A 2019, 52, 425601. [Google Scholar] [CrossRef]
  8. Sene, N. SIR epidemic model with Mittag-Leffler fractional derivative. Chaos Solitons Fractals 2020, 137, 109833. [Google Scholar] [CrossRef]
  9. Simon, M. SIR epidemics with stochastic infectious periods. Stoch. Proc. Appl. 2020, 130, 4252–4274. [Google Scholar] [CrossRef]
  10. Tian, C.R.; Zhang, Q.Y.; Zhang, L. Global stability in a networked SIR epidemic model. Appl. Math. Lett. 2020, 107, 106444. [Google Scholar] [CrossRef]
  11. El Koufi, A.; Adnani, J.; Bennar, A.; Yousfi, N. Analysis of a stochastic SIR model with vaccination and nonlinear incidence rate. Int. J. Diff. Equ. 2019, 2019, 9275051. [Google Scholar] [CrossRef] [Green Version]
  12. Houy, N. Are better vaccines really better? The case of a simple stochastic epidemic SIR model. Econ. Bull. 2013, 33, 207–216. [Google Scholar]
  13. Jornet-Sanz, M.; Corberan-Vallet, A.; Santonja, F.J.; Villanueva, R.J. A bayesian stochastic SIRS model with a vaccination strategy for the analysis of respiratory syncytial virus. Stat. Oper. Res. Trans. 2017, 41, 159–175. [Google Scholar]
  14. Li, X.N.; Zhang, Q.M. Time to extinction and stationary distribution of a stochastic susceptible-infected-recovered-susceptible model with vaccination under markov switching. Math. Popul. Stud. 2020, 27, 259–274. [Google Scholar] [CrossRef]
  15. Liu, Q.; Jiang, D.Q. The threshold of a stochastic delayed SIR epidemic model with vaccination. Phys. A 2016, 461, 140–147. [Google Scholar] [CrossRef]
  16. Liu, Q.; Jiang, D.Q.; Shi, N.Z.; Hayat, T. Dynamics of a stochastic delayed SIR epidemic model with vaccination and double diseases driven by levy jumps. Phys. A 2018, 492, 2010–2018. [Google Scholar] [CrossRef]
  17. Miao, A.Q.; Zhang, J.; Zhang, T.Q.; Pradeep, B. Threshold dynamics of a stochastic SIR model with vertical transmission and vaccination. Comput. Math. Meth. Med. 2017, 2017, 4820183. [Google Scholar] [CrossRef] [Green Version]
  18. Nguyen, C.; Carlson, J.M. Optimizing real-time vaccine allocation in a stochastic SIR model. PLoS ONE 2016, 11, 0152950. [Google Scholar] [CrossRef]
  19. Wang, F.Y.; Wang, X.Y.; Zhang, S.W.; Ding, C.M. On pulse vaccine strategy in a periodic stochastic SIR epidemic model. Chaos Solitons Fractals 2014, 66, 127–135. [Google Scholar] [CrossRef]
  20. Wang, L.; Teng, Z.D.; Tang, T.T.; Li, Z.M. Threshold dynamics in stochastic SIRS epidemic models with nonlinear incidence and vaccination. Comput. Math. Meth. Med. 2017, 2017, 7294761. [Google Scholar] [CrossRef] [Green Version]
  21. Witbooi, P.J. Stability of a stochastic model of an SIR epidemic with vaccination. Acta Biotheor. 2017, 65, 151–165. [Google Scholar] [CrossRef] [PubMed]
  22. Xu, C.Y.; Li, X.Y. The threshold of a stochastic delayed SIRS epidemic model with temporary immunity and vaccination. Chaos Solitons Fractals 2018, 111, 227–234. [Google Scholar] [CrossRef]
  23. Zhang, Y.; Li, Y.; Zhang, Q.L.; Li, A.H. Behavior of a stochastic SIR epidemic model with saturated incidence and vaccination rules. Phys. A 2018, 501, 178–187. [Google Scholar] [CrossRef]
  24. Zhao, X.; He, X.; Feng, T.; Qiu, Z.P. A stochastic switched SIRS epidemic model with nonlinear incidence and vaccination: Stationary distribution and extinction. Int. J. Biomath. 2020, 13, 2050020. [Google Scholar] [CrossRef]
  25. Colombo, R.M.; Garavello, M. Optimizing vaccination strategies in an age structured SIR model. Math. Biosci. Eng. 2020, 17, 1074–1089. [Google Scholar] [CrossRef]
  26. Cui, Q.Q.; Xu, J.B.; Zhang, Q.; Wang, K. An nsfd scheme for sir epidemic models of childhood diseases with constant vaccination strategy. Adv. Differ. Equ. 2014, 2014, 172. [Google Scholar] [CrossRef] [Green Version]
  27. D’Onofrio, A. Pulse vaccination strategy in the SIR epidemic model: Global asymptotic stable eradication in presence of vaccine failures. Math. Comput. Model. 2002, 36, 473–489. [Google Scholar] [CrossRef]
  28. d’Onofrio, A. On pulse vaccination strategy in the SIR epidemic model with vertical transmission. Appl. Math. Lett. 2005, 18, 729–732. [Google Scholar] [CrossRef] [Green Version]
  29. Gao, S.J.; Ouyang, H.S.; Nieto, J.J. Mixed vaccination strategy in SIRS epidemic model with seasonal variability on infection. Int. J. Biomath. 2011, 4, 473–491. [Google Scholar] [CrossRef]
  30. Gao, S.J.; Xie, D.H.; Chen, L.S. Pulse vaccination strategy in a delayed SIR epidemic model with vertical transmission. Discr. Contin. Dyn. Syst. B 2007, 7, 77–86. [Google Scholar] [CrossRef]
  31. Kabir, K.M.A.; Tanimoto, J. Vaccination strategies in a two-layer sir/v-ua epidemic model with costly information and buzz effect. Commun. Nonlin. Sci. Numer. Simul. 2019, 76, 92–108. [Google Scholar] [CrossRef]
  32. Li, J.Q.; Yang, Y.L. SIR-SVS epidemic models with continuous and impulsive vaccination strategies. J. Theor. Biol. 2011, 280, 108–116. [Google Scholar] [CrossRef]
  33. Liu, H.L.; Yu, J.Y.; Zhu, G.T. Global stability of an age-structured SIR epidemic model with pulse vaccination strategy. J. Syst. Sci. Complex. 2012, 25, 417–429. [Google Scholar] [CrossRef]
  34. Liu, L.; Luo, X.F.; Chang, L.L. Vaccination strategies of an SIR pair approximation model with demographics on complex networks. Chaos Solitons Fractals 2017, 104, 282–290. [Google Scholar] [CrossRef]
  35. Makinde, O.D. Adomian decomposition approach to a SIR epidemic model with constant vaccination strategy. Appl. Math. Comput. 2007, 184, 842–848. [Google Scholar] [CrossRef]
  36. Meng, X.Z.; Chen, L.S. The dynamics of a new SIR epidemic model concerning pulse vaccination strategy. Appl. Math. Comput. 2008, 197, 582–597. [Google Scholar] [CrossRef]
  37. Moneim, I.A.; Greenhalgh, D. Threshold and stability results for an SIRS epidemic model with a general periodic vaccination strategy. J. Biol. Syst. 2005, 13, 131–150. [Google Scholar] [CrossRef]
  38. Mu, X.J.; Zhang, Q.M.; Rong, L.B. Optimal vaccination strategy for an SIRS model with imprecise parameters and levy noise. J. Franklin Inst. Eng. Appl. Math. 2019, 356, 11385–11413. [Google Scholar] [CrossRef]
  39. Mungkasi, S. Variational iteration and successive approximation methods for a SIR epidemic model with constant vaccination strategy. Appl. Math. Model. 2021, 90, 1–10. [Google Scholar] [CrossRef]
  40. Pei, Y.Z.; Liu, S.Y.; Chen, L.S.; Wang, C.H. Two different vaccination strategies in an SIR epidemic model with saturated infectious force. Int. J. Biomath. 2008, 1, 147–160. [Google Scholar] [CrossRef]
  41. Rashidinia, J.; Sajjadian, M.; Duarte, J.; Januario, C.; Martins, N. On the dynamical complexity of a seasonally forced discrete SIR epidemic model with a constant vaccination strategy. Complexity 2018, 2018, 7191487. [Google Scholar] [CrossRef]
  42. Terry, A.J. PULSE vaccination strategies in a metapopulation SIR model. Math. Biosci. Eng. 2010, 7, 455–477. [Google Scholar] [PubMed]
  43. Zhou, L.H.; Wang, Y.; Xiao, Y.Y.; Li, M.Y. Global dynamics of a discrete age-structured SIR epidemic model with applications to measles vaccination strategies. Math. Biosci. 2019, 308, 27–37. [Google Scholar] [CrossRef] [PubMed]
  44. Karatayev, V.A.; Anand, M.; Bauch, C.T. Local lockdowns outperform global lockdown on the far side of the COVID-19 epidemic curve. Proc. Natl. Acad. Sci. USA 2020, 117, 24575–24580. [Google Scholar] [CrossRef]
  45. Wang, X.; Jia, D.; Gao, S.; Xia, C.; Li, X.; Wang, Z. Vaccination behavior by coupling the epidemic spreading with the human decision under the game theory. Appl. Math. Comput. 2020, 380, 125232. [Google Scholar] [CrossRef]
  46. Priesemann, V.; Balling, R.; Brinkmann, M.M.; Ciesek, S.; Czypionka, T.; Eckerle, I.; Giordano, G.; Hanson, C.; Hel, Z.; Hotulainen, P.; et al. An action plan for pan-European defence against new SARS-CoV-2 variants. Lancet 2021, 397, 469–470. [Google Scholar] [CrossRef]
  47. Priesemann, V.; Brinkmann, M.M.; Ciesek, S.; Cuschieri, S.; Czypionka, T.; Giordano, G.; Gurdasani, D.; Hanson, C.; Hens, N.; Iftekhar, E.; et al. Calling for pan-European commitment for rapid and sustained reduction in SARS-CoV-2 infections. Lancet 2021, 397, 92–93. [Google Scholar] [CrossRef]
  48. Zhou, Y.G.; Yang, K.; Zhou, K.; Liang, Y.T. Optimal vaccination policies for an SIR model with limited resources. Acta Biotheor. 2014, 62, 171–181. [Google Scholar] [CrossRef]
  49. Abouelkheir, I.; El Kihal, F.; Rachik, M.; Elmouki, I. Optimal impulse vaccination approach for an SIR control model with short-term immunity. Mathematics 2019, 7, 420. [Google Scholar] [CrossRef] [Green Version]
  50. Church, K.E.M.; Liu, X.Z. Analysis of a SIR model with pulse vaccination and temporary immunity: Stability, bifurcation and a cylindrical attractor. Nonlinear Anal. Real World Appl. 2019, 50, 240–266. [Google Scholar] [CrossRef]
  51. Gao, S.J.; Teng, Z.D.; Xie, D.H. Analysis of a delayed SIR epidemic model with pulse vaccination. Chaos Solitons Fractals 2009, 40, 1004–1011. [Google Scholar] [CrossRef]
  52. Gao, S.J.; Zhidong, Z.D.; Nieto, J.J.; Torres, A. Analysis of an SIR epidemic model with pulse vaccination and distributed time delay. J. Biomed. Biotechn. 2007, 2007, 64870. [Google Scholar] [CrossRef]
  53. He, Z.L.; Nie, L.F. The effect of pulse vaccination and treatment on SIR epidemic model with media impact. Discr. Dyn. Nat. Soc. 2015, 2015, 532494. [Google Scholar] [CrossRef] [Green Version]
  54. Jiang, G.R.; Yang, Q.G. Bifurcation analysis in an SIR epidemic model with birth pulse and pulse vaccination. Appl. Math. Comput. 2009, 215, 1035–1046. [Google Scholar] [CrossRef]
  55. Liu, X.S.; Dai, B.X. Qualitative and bifurcation analysis of an SIR epidemic model with saturated treatment function and nonlinear pulse vaccination. Discr. Dyn. Nat. Soc. 2016, 2016, 9146481. [Google Scholar] [CrossRef]
  56. Liu, X.S.; Dai, B.X. Flip bifurcations of an SIR epidemic model with birth pulse and pulse vaccination. Appl. Math. Model. 2017, 43, 579–591. [Google Scholar] [CrossRef]
  57. Lu, Z.H.; Chi, X.B.; Chen, L.S. The effect of constant and pulse vaccination on SIR epidemic model with horizontal and vertical transmission. Math. Comput. Model. 2002, 36, 1039–1057. [Google Scholar] [CrossRef]
  58. Meng, X.Z.; Chen, L.S. Global dynamical behaviors for an SIR epidemic model with time delay and pulse vaccination. Taiwan. J. Math. 2008, 12, 1107–1122. [Google Scholar] [CrossRef]
  59. Meng, X.Z.; Chen, L.S.; Wu, B. A delay SIR epidemic model with pulse vaccination and incubation times. Nonlin. Anal. Real World Appl. 2010, 11, 88–98. [Google Scholar] [CrossRef]
  60. Nie, L.F.; Teng, Z.D.; Torres, A. Dynamic analysis of an SIR epidemic model with state dependent pulse vaccination. Nonlin. Anal. Real World Appl. 2012, 13, 1621–1629. [Google Scholar] [CrossRef]
  61. Pang, G.P.; Chen, L.S. A delayed SIRS epidemic model with pulse vaccination. Chaos Solitons Fractals 2007, 34, 1629–1635. [Google Scholar] [CrossRef]
  62. Qin, W.J.; Tang, S.Y.; Cheke, R.A. Nonlinear pulse vaccination in an SIR epidemic model with resource limitation. Abstr. Appl. Anal. 2013, 2013, 670263. [Google Scholar] [CrossRef]
  63. Sekiguchi, M.; Ishiwata, E. Dynamics of a discretized SIR epidemic model with pulse vaccination and time delay. J. Comput. Appl. Math. 2011, 236, 997–1008. [Google Scholar] [CrossRef] [Green Version]
  64. Stone, L.; Shulgin, B.; Agur, Z. Theoretical examination of the pulse vaccination policy in the SIR epidemic model. Math. Comput. Model. 2000, 31, 207–215. [Google Scholar] [CrossRef]
  65. Wang, L. Existence of periodic solutions of seasonally forced SIR models with impulse vaccination. Taiwan. J. Math. 2015, 19, 1713–1729. [Google Scholar] [CrossRef]
  66. Zhang, X.B.; Huo, H.F.; Sun, X.K.; Fu, Q. The differential susceptibility SIR epidemic model with time delay and pulse vaccination. J. Appl. Math. Comput. 2010, 34, 287–298. [Google Scholar] [CrossRef]
  67. Zhang, X.B.; Huo, H.F.; Sun, X.K.; Fu, Q.A. The differential susceptibility SIR epidemic model with stage structure and pulse vaccination. Nonlin. Anal. Real World Appl. 2010, 11, 2634–2646. [Google Scholar] [CrossRef]
  68. Zhang, X.B.; Huo, H.F.; Xiang, H.; Meng, X.Y. An SIRS epidemic model with pulse vaccination and non-monotonic incidence rate. Nonlin. Anal. Hybrid Syst. 2013, 8, 13–21. [Google Scholar] [CrossRef]
  69. Zhao, W.C.; Li, J.; Meng, X.Z. Dynamical analysis of SIR epidemic model with nonlinear pulse vaccination and lifelong immunity. Discret. Dyn. Nat. Soc. 2015, 2015, 848623. [Google Scholar] [CrossRef]
  70. Zhao, Z.; Pang, L.Y.; Chen, Y. Nonsynchronous bifurcation of SIRS epidemic model with birth pulse and pulse vaccination. Nonlin. Dyn. 2015, 79, 2371–2383. [Google Scholar] [CrossRef]
  71. Zhou, A.R.; Sattayatham, P.; Jiao, J.J. Dynamics of an SIR epidemic model with stage structure and pulse vaccination. Adv. Diff. Equ. 2016, 2016, 140. [Google Scholar] [CrossRef] [Green Version]
  72. d’Onofrio, A.; Manfredi, P.; Salinelli, E. Bifurcation thresholds in an SIR model with information-dependent vaccination. Math. Model. Nat. Phenom. 2007, 2, 26–43. [Google Scholar] [CrossRef]
  73. Gumus, O.A.K.; Selvam, A.G.M.; Vianny, D.A. Bifucaction and stability analysis of a discrete time SIR epidemic model with vaccination. Int. J. Anal. Appl. 2019, 17, 809–820. [Google Scholar]
  74. Rostamy, D.; Mottaghi, E. Forward and backward bifurcation in a fractional-order SIR epidemic model with vaccination. Iran. J. Sci. Technol. Trans. A 2018, 42, 663–671. [Google Scholar] [CrossRef]
  75. Zhang, Q.Q.; Tang, B.; Tang, S.Y. Vaccination threshold size and backward bifurcation of SIR model with state-dependent pulse control. J. Theor. Biol. 2018, 455, 75–85. [Google Scholar] [CrossRef]
  76. Elazzouzi, A.; Alaoui, A.L.; Tilioua, M.; Tridane, A. Global stability analysis for a generalized delayed SIR model with vaccination and treatment. Adv. Diff. Equ. 2019, 2019, 532. [Google Scholar] [CrossRef]
  77. Laarabi, H.; Abta, A.; Hattaf, K. Optimal control of a delayed SIRS epidemic model with vaccination and treatment. Acta Biotheor. 2015, 63, 87–97. [Google Scholar] [CrossRef]
  78. Liu, Q.; Jiang, D.Q.; Hayat, T.; Ahmad, B. Analysis of a delayed vaccinated SIR epidemic model with temporary immunity and levy jumps. Nonlin. Anal. Hybrid Syst. 2018, 27, 29–43. [Google Scholar] [CrossRef]
  79. Tian, X.H. Stability analysis of a delayed SIRS epidemic model with vaccination and nonlinear incidence. Int. J. Biomath. 2012, 5, 1250050. [Google Scholar] [CrossRef]
  80. Bakare, E.A. On the optimal control of vaccination and treatments for an SIR-epidemic model with infected immigrants. Int. J. Ecol. Econ. Statist. 2016, 37, 82–104. [Google Scholar]
  81. Chapman, J.D.; Evans, N.D. The structural identifiability of susceptible-infective-recovered type epidemic models with incomplete immunity and birth targeted vaccination. Biomed. Signal Proc. Control 2009, 4, 278–284. [Google Scholar] [CrossRef]
  82. Guo, H.J.; Chen, L.S.; Song, X.Y. Dynamical properties of a kind of SIR model with constant vaccination rate and impulsive state feedback control. Int. J. Biomath. 2017, 10, 17500930. [Google Scholar] [CrossRef]
  83. Kar, T.K.; Batabyal, A. Stability analysis and optimal control of an SIR epidemic model with vaccination. Biosystems 2011, 104, 127–135. [Google Scholar] [CrossRef]
  84. Ledzewicz, U.; Schattler, H. ON optimal singular controls for a general sir-model with vaccination and treatment. Discr. Contin. Dyn. Syst. 2011, 31, 981–990. [Google Scholar]
  85. Rao, F.; Mandal, P.S.; Kang, Y. Complicated endemics of an SIRS model with a generalized incidence under preventive vaccination and treatment controls. Appl. Math. Model. 2019, 67, 38–61. [Google Scholar] [CrossRef]
  86. Rao, X.B.; Zhao, X.P.; Chu, Y.D.; Zhang, J.G.; Gao, J.S. The analysis of mode-locking topology in an SIR epidemic dynamics model with impulsive vaccination control: Infinite cascade of stern-brocot sum trees. Chaos Solitons Fractals 2020, 139, 110031. [Google Scholar] [CrossRef]
  87. Zeng, G.Z.; Chen, L.S.; Sun, L.H. Complexity of an SIR epidemic dynamics model with impulsive vaccination control. Chaos Solitons Fractals 2005, 26, 495–505. [Google Scholar] [CrossRef]
  88. Elbasha, E.H.; Podder, C.N.; Gumel, A.B. Analyzing the dynamics of an SIRS vaccination model with waning natural and vaccine-induced immunity. Nonlin. Anal. Real World Appl. 2011, 12, 2692–2705. [Google Scholar] [CrossRef]
  89. Hui, J.; Chen, L.S. Impulsive vaccination of SIR epidemic models with nonlinear incidence rates. Discr. Contin. Dyn. Syst. B 2004, 4, 595–605. [Google Scholar] [CrossRef]
  90. Khader, M.M.; Adel, M. Numerical treatment of the fractional modeling on susceptible-infected-recovered equations with a constant vaccination rate by using gem. Int. J. Nonlin. Sci. Numer. Simul. 2019, 20, 69–75. [Google Scholar] [CrossRef]
  91. Li, B.L.; Qin, C.Y.; Wang, X. Analysis of an SIRS epidemic model with nonlinear incidence and vaccination. Commun. Math. Biol. Neurosci. 2020, 2020, 4262. [Google Scholar]
  92. Sun, C.J.; Yang, W. Global results for an SIRS model with vaccination and isolation. Nonlin. Anal. Real World Appl. 2010, 11, 4223–4237. [Google Scholar] [CrossRef]
  93. Lahrouz, A.; Omari, L.; Kiouach, D.; Belmaati, A. Complete global stability for an SIRS epidemic model with generalized non-linear incidence and vaccination. Appl. Math. Comput. 2012, 218, 6519–6525. [Google Scholar] [CrossRef]
  94. Wang, X.; Gao, S.; Zhu, P.; Wang, J. Roles of different update strategies in the vaccination behavior on two-layered networks. Phys. Lett. A 2020, 384, 126224. [Google Scholar] [CrossRef]
  95. Assadouq, A.; El Mahjour, H.; Settati, A. Qualitative behavior of a SIRS epidemic model with vaccination on heterogeneous networks. Ital. J. Pure Appl. Math. 2020, 43, 958–974. [Google Scholar]
  96. Le Chang, S.; Piraveenan, M.; Prokopenko, M. The effects of imitation dynamics on vaccination behaviours in sir-network model. Int. J. Environ. Res. Public Health 2019, 16, 16142477. [Google Scholar] [CrossRef] [Green Version]
  97. Auchincloss, A.H.; Roux, A.V.D. A new tool for epidemiology: The usefulness of dynamic-agent models in understanding place effects on health. Am. J. Epidemiol. 2008, 168, 1–8. [Google Scholar] [CrossRef] [Green Version]
  98. Ajelli, M.; Goncalves, B.; Balcan, D.; Colizza, V.; Hu, H.; Ramasco, J.J.; Merler, S.; Vespignani, A. Comparing large-scale computational approaches to epidemic modeling: Agent-based versus structured metapopulation models. BMC Infect. Dis. 2010, 10, 190. [Google Scholar] [CrossRef] [Green Version]
  99. Schüttler, J.; Schlickeiser, R.; Schlickeiser, F.; Kröger, M. Covid-19 predictions using a Gauss model, based on data from April 2. Physics 2020, 2, 197–212. [Google Scholar] [CrossRef]
  100. Yildirim, A.; Cherruault, Y. Analytical approximate solution of a SIR epidemic model with constant vaccination strategy by homotopy perturbation method. Kybernetes 2009, 38, 1566–1575. [Google Scholar] [CrossRef]
  101. Heng, K.; Althaus, C.L. The approximately universal shapes of epidemic curves in the susceptible-exposed-infectious-recovered (SEIR) model. Sci. Rep. 2020, 10, 19365. [Google Scholar] [CrossRef]
  102. Barlow, N.S.; Weinstein, S.J. Accurate closed-form solution of the SIR epidemic model. Phys. D 2020, 408, 132540. [Google Scholar] [CrossRef]
  103. Bidari, S.; Chen, X.Y.; Peters, D.; Pittman, D.; Simon, P.L. Solvability of implicit final size equations for SIR epidemic models. Math. Biosci. 2016, 282, 181–190. [Google Scholar] [CrossRef] [Green Version]
  104. Carvalho, A.M.; Goncalves, S. An analytical solution for the Kermack-McKendrick model. Phys. A 2021, 566, 125659. [Google Scholar] [CrossRef]
  105. Guerrero, F.; Santonja, F.J.; Villanueva, R.J. Solving a model for the evolution of smoking habit in Spain with homotopy analysis method. Nonlinear Anal. Real World Appl. 2013, 14, 549–558. [Google Scholar] [CrossRef]
  106. Khan, H.; Mohapatra, R.N.; Vajravelu, K.; Liao, S.J. The explicit series solution of SIR and SIS epidemic models. Appl. Math. Comput. 2009, 215, 653–669. [Google Scholar] [CrossRef]
  107. Liu, J.L.; Peng, B.Y.; Zhang, T.L. Effect of discretization on dynamical behavior of SEIR and SIR models with nonlinear incidence. Appl. Math. Lett. 2015, 39, 60–66. [Google Scholar] [CrossRef]
  108. Van Mieghem, P. Approximate formula and bounds for the time-varying susceptible-infected-susceptible prevalence in networks. Phys. Rev. E 2016, 93, 052312. [Google Scholar] [CrossRef]
  109. Kröger, M.; Schlickeiser, R. Analytical solution of the SIR-model for the temporal evolution of epidemics. Part A: Time-independent reproduction factor. J. Phys. A Math. Theor. 2020, 53, 505601. [Google Scholar] [CrossRef]
  110. Schlickeiser, R.; Kröger, M. Analytical solution of the SIR-model for the temporal evolution of epidemics. Part B: Semi-time case. J. Phys. A Math. Theor. 2021, 54, 175601. [Google Scholar] [CrossRef]
  111. Turkyilmazoglu, M. Explicit formulae for the peak time of an epidemic from the SIR model. Phys. D 2021, 422, 132902. [Google Scholar] [CrossRef] [PubMed]
  112. Estrada, E. Covid-19 and Sars-Cov-2. Modeling the present, looking at the future. Phys. Rep. 2020, 869, 1. [Google Scholar] [CrossRef] [PubMed]
  113. Kröger, M.; Schlickeiser, R. Forecast for the second Covid-19 wave based on the improved SIR-model with a constant ratio of recovery to infection rate. Preprints 2021, 2021010449. [Google Scholar] [CrossRef]
  114. Morton, R.; Wickwire, K.H. On the optimal control of a deterministic epidemic. Adv. Appl. Probab. 1974, 6, 622–635. [Google Scholar] [CrossRef]
  115. Anderson, R.M.; May, R.M. Infectious Diseases of Humans: Dynamics and Control; Oxford Science Publications: Oxford, UK, 1991. [Google Scholar]
  116. Behncke, H. Optimal control of deterministic epidemics. Optim. Control Appl. Methods 2000, 21, 269–285. [Google Scholar] [CrossRef]
  117. Hansen, E.; Day, T. Optimal control of epidemics with limited resources. J. Math. Biol. 2011, 62, 423–451. [Google Scholar] [CrossRef] [Green Version]
  118. Grauer, J.; Löwen, H.; Liebchen, B. Strategic spatiotemporal vaccine distribution increases the survival rate in an infectious disease like Covid-19. Sci. Rep. 2021, 10, 21594. [Google Scholar] [CrossRef]
  119. Grundel, S.; Heyder, S.; Hotz, T.; Ritschel, T.K.S.; Sauerteig, P.; Worthmann, K. How to coordinate vaccination and social distancing to mitigate SARS-CoV-2 outbreaks. medRxiv 2020, 2020, 20248707. [Google Scholar] [CrossRef]
  120. Duclos, T.; Reichert, T. The missing Link: A closed form solution to the Kermack and McKendrick epidemic model equations. medRxiv 2021, 2021, 21252781. [Google Scholar] [CrossRef]
  121. Abramowitz, M.; Stegun, I.A. Handbook of Mathematical Functions; National Bureau of Standards: Washington, DC, USA, 1972.
  122. Kröger, M.; Schlickeiser, R. Gaussian doubling times and reproduction factors of the COVID-19 pandemic disease. Front. Phys. 2020, 8, 276. [Google Scholar] [CrossRef]
  123. Data Repository. 2021. Available online: https://github.com/owid/covid-19-data/blob/master/public/data/vaccinations/vaccinations.csv (accessed on 20 May 2021).
  124. Data Repository. 2021. Available online: https://www.complexfluids.ethz.ch/cgi-bin/covid19-waveII (accessed on 20 May 2021).
Figure 1. Three time-dependent rates a ( t ) , μ ( t ) , and v ( t ) entering the SIRV equations for the four compartments of susceptible, S, infectious, I, recovered, R, and vaccinated, V, population fractions. Upon introducing reduced time τ , the model is characterized by the assumed constant ratios k = μ ( t ) / a ( t ) and b = v ( t ) / a ( t ) .
Figure 1. Three time-dependent rates a ( t ) , μ ( t ) , and v ( t ) entering the SIRV equations for the four compartments of susceptible, S, infectious, I, recovered, R, and vaccinated, V, population fractions. Upon introducing reduced time τ , the model is characterized by the assumed constant ratios k = μ ( t ) / a ( t ) and b = v ( t ) / a ( t ) .
Physics 03 00028 g001
Figure 2. Exact ψ m versus k and b in the lower right triangle, and approximate ψ m using Equation (74) above the diagonal (mirrored, to allow for a simple comparison with the exact ψ m ). All analytic results for the SIRV functions in terms of reduced time τ are basically exact for those k and b for which ψ m is well described by its approximant. The white space in the lower left corner is the regime of b < b c , where the SIRV model results are captured by linearly interpolating between the SIR model and the SIRV model evaluated at the critical b = b c .
Figure 2. Exact ψ m versus k and b in the lower right triangle, and approximate ψ m using Equation (74) above the diagonal (mirrored, to allow for a simple comparison with the exact ψ m ). All analytic results for the SIRV functions in terms of reduced time τ are basically exact for those k and b for which ψ m is well described by its approximant. The white space in the lower left corner is the regime of b < b c , where the SIRV model results are captured by linearly interpolating between the SIR model and the SIRV model evaluated at the critical b = b c .
Physics 03 00028 g002
Figure 3. Comparison of the approximant (80) for ψ ( τ ) / ψ 0 with the exact solutions (59) and (61) (black curves) for three different α values at η = 10 6 and a relatively low b = 0.02 . For the approximant (green), ψ m and τ m are given by Equations (71) and (74), respectively. For larger b the performance of the approximant is even better.
Figure 3. Comparison of the approximant (80) for ψ ( τ ) / ψ 0 with the exact solutions (59) and (61) (black curves) for three different α values at η = 10 6 and a relatively low b = 0.02 . For the approximant (green), ψ m and τ m are given by Equations (71) and (74), respectively. For larger b the performance of the approximant is even better.
Physics 03 00028 g003
Figure 4. The critical b c versus k and η . The coloring scheme uses the decadic logarithm of b c , and the vertical axis is also logarithmic. Only a relevant range of k values is shown.
Figure 4. The critical b c versus k and η . The coloring scheme uses the decadic logarithm of b c , and the vertical axis is also logarithmic. Only a relevant range of k values is shown.
Physics 03 00028 g004
Figure 5. Reduced peak time, τ max j , of the newly infected population fraction, j ( τ ) , versus reduced vaccination rate, b, for various k at η = 10 6 (double-logarithmic plot). The exact numerical solution (solid black) is compared with the approximant (98) within the regime of b > b c (thin colored), and by the corresponding linear interpolant (thick colored) for the remaining regime of very small b < b c . The limiting value τ max j ( b 0 ) exactly coincides with the τ max SIR ( k ) of the SIR model (see Appendix G.2).
Figure 5. Reduced peak time, τ max j , of the newly infected population fraction, j ( τ ) , versus reduced vaccination rate, b, for various k at η = 10 6 (double-logarithmic plot). The exact numerical solution (solid black) is compared with the approximant (98) within the regime of b > b c (thin colored), and by the corresponding linear interpolant (thick colored) for the remaining regime of very small b < b c . The limiting value τ max j ( b 0 ) exactly coincides with the τ max SIR ( k ) of the SIR model (see Appendix G.2).
Physics 03 00028 g005
Figure 6. Peak value of the newly infected population fraction j max versus reduced vaccination rate, b, for various k at η = 10 6 (double-logarithmic plot). The exact numerical solution (solid black) is compared with the approximant (99) within the regime of b > b c (thin colored), and by the interpolant (101) (thick colored) for the remaining regime of very small b < b c . The limiting value j max ( b 0 ) exactly coincides with the j max SIR ( k ) of the SIR model (Appendix G.2).
Figure 6. Peak value of the newly infected population fraction j max versus reduced vaccination rate, b, for various k at η = 10 6 (double-logarithmic plot). The exact numerical solution (solid black) is compared with the approximant (99) within the regime of b > b c (thin colored), and by the interpolant (101) (thick colored) for the remaining regime of very small b < b c . The limiting value j max ( b 0 ) exactly coincides with the j max SIR ( k ) of the SIR model (Appendix G.2).
Physics 03 00028 g006
Figure 7. Final (infinite time) fraction of infected persons, J , versus reduced vaccination rate, b, for various k at η = 10 6 (double-logarithmic plot). The exact numerical solution (solid black) is compared with the approximant (102) within the regime of b > b c (thin colored), and by the interpolant (103) (thick colored) for the remaining regime of very small b < b c . The limiting value J ( b 0 ) exactly coincides with the J SIR ( k ) of the SIR model (see Appendix G.2).
Figure 7. Final (infinite time) fraction of infected persons, J , versus reduced vaccination rate, b, for various k at η = 10 6 (double-logarithmic plot). The exact numerical solution (solid black) is compared with the approximant (102) within the regime of b > b c (thin colored), and by the interpolant (103) (thick colored) for the remaining regime of very small b < b c . The limiting value J ( b 0 ) exactly coincides with the J SIR ( k ) of the SIR model (see Appendix G.2).
Physics 03 00028 g007
Figure 8. Differential rate, j ( τ ) , of infected population fraction versus reduced time, τ , for three different k { 0.85 , 0.90 , 0.95 } and various reduced vaccination rates, b / b c . Here, η = 10 6 is used. The panels (ac) show the regime b > b c , while (df) show results for b < b c including b = 0 (SIR model). To rate the effect of the parameters, all three plots of each row are shown on identical scales. While the peak time increases with increasing k and decreasing b, the peak height dramatically decreases with increasing b. The critical b c depends on k and η , see Equation (92) and Figure 4. The area under the curves is the total cumulative fraction J of infected persons. The black lines are the exact numerical results, the green lines are the analytical approximant, given in Section 4.7.
Figure 8. Differential rate, j ( τ ) , of infected population fraction versus reduced time, τ , for three different k { 0.85 , 0.90 , 0.95 } and various reduced vaccination rates, b / b c . Here, η = 10 6 is used. The panels (ac) show the regime b > b c , while (df) show results for b < b c including b = 0 (SIR model). To rate the effect of the parameters, all three plots of each row are shown on identical scales. While the peak time increases with increasing k and decreasing b, the peak height dramatically decreases with increasing b. The critical b c depends on k and η , see Equation (92) and Figure 4. The area under the curves is the total cumulative fraction J of infected persons. The black lines are the exact numerical results, the green lines are the analytical approximant, given in Section 4.7.
Physics 03 00028 g008
Figure 9. Suitable normalized SIRV quantities S, I / I max , R / R , V, and J / J versus τ for the four different reduced vaccination rates at k = 0.9 and η = 10 6 : (a,e) b = 0 (SIR model), (b,f) b = 0.1 b c , (c,g) b = b c , and (d,h) b = 10 b c . While the upper row presents the data in a linear scale, the bottom row shows the same but in a semilogarithmic fashion to appreciate the two well separated time scales.
Figure 9. Suitable normalized SIRV quantities S, I / I max , R / R , V, and J / J versus τ for the four different reduced vaccination rates at k = 0.9 and η = 10 6 : (a,e) b = 0 (SIR model), (b,f) b = 0.1 b c , (c,g) b = b c , and (d,h) b = 10 b c . While the upper row presents the data in a linear scale, the bottom row shows the same but in a semilogarithmic fashion to appreciate the two well separated time scales.
Physics 03 00028 g009
Figure 10. Same as Figure 4 but using another colormap, where countries have been added. 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 Table 1, and the brightness of a filled circle corresponds to the reduced vaccination rate b, also taken from Table 1.
Figure 10. Same as Figure 4 but using another colormap, where countries have been added. 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 Table 1, and the brightness of a filled circle corresponds to the reduced vaccination rate b, also taken from Table 1.
Physics 03 00028 g010
Table 1. Analysis using data from 18 March 2021. For η , k and a the current values for the second wave, that started at t 0 II are used, all from the online resource [124]. The starting time, t V , of the vaccination program and the mean daily fraction v of vaccinated population since then are retrieved from Ref. [123] assuming that each person has to be vaccinated twice and that the vaccination is effective two weeks after the second shot. The remaining quantities are derived from Equations (9) and (36), i.e., via μ = a k , b = v / a , α = k b , and Δ = 1 2 η k b is positive if the outburst condition (12) is fulfilled. b c is calculated via Equation (92). Furthermore, included are the infected population fraction at various times: (i) J I at the end of the first wave, J ( t V ) at the onset of vaccinations, (iii) J b = 0 assuming no vaccinations, (iv) J b = b assuming ongoing vaccination at the present rate, (v) J b = 2 b 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, and J ˙ 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 our Supplementary Material.
Table 1. Analysis using data from 18 March 2021. For η , k and a the current values for the second wave, that started at t 0 II are used, all from the online resource [124]. The starting time, t V , of the vaccination program and the mean daily fraction v of vaccinated population since then are retrieved from Ref. [123] assuming that each person has to be vaccinated twice and that the vaccination is effective two weeks after the second shot. The remaining quantities are derived from Equations (9) and (36), i.e., via μ = a k , b = v / a , α = k b , and Δ = 1 2 η k b is positive if the outburst condition (12) is fulfilled. b c is calculated via Equation (92). Furthermore, included are the infected population fraction at various times: (i) J I at the end of the first wave, J ( t V ) at the onset of vaccinations, (iii) J b = 0 assuming no vaccinations, (iv) J b = b assuming ongoing vaccination at the present rate, (v) J b = 2 b 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, and J ˙ 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 our Supplementary Material.
Countryka μ b b / b c α Δ t 0 II t V J I J ( t V ) J b = 0 J b = b J b = 2 b J ˙ max t 99 %
α 3 Code [d 1 ][d 1 ]
ARG0.9120.1250.1140.00220.0390.9100.0720-08-1720-12-280.0690.2280.2820.2770.27412921-10-17
AUT0.9050.5200.4710.00130.8410.9040.0920-07-2920-12-260.0080.1280.1910.1820.17721921-07-29
BEL0.8960.5510.4940.00110.2710.8950.1020-09-0620-12-270.1620.3130.3320.3300.32924421-05-26
BRA0.7900.0460.0370.00990.2020.7800.1920-07-0521-01-140.1450.2470.4890.4320.3988622-09-29
CAN0.9621.0180.9800.00040.6920.9620.0420-09-2620-12-130.0490.0690.1210.1060.0996821-06-05
CHE0.8940.4580.4090.00231.2370.8910.1020-07-2221-01-220.0440.2160.2400.2360.23423121-09-05
DEU0.9150.5590.5110.00121.0770.9130.0820-08-1620-12-260.0190.0690.1820.1580.14318021-08-06
ESP0.8760.1750.1530.00461.1200.8710.1220-04-2921-01-020.0920.2070.3090.2830.27011522-01-28
FIN0.9973.8583.8480.00022.1810.9970.0020-12-2120-12-300.0070.0080.0190.0120.0111721-02-21
FRA0.8860.2280.2020.00270.9730.8830.1120-05-1120-12-260.0820.1830.2840.2630.25212721-12-25
GBR0.8670.3890.3370.00531.7410.8620.1320-09-0920-12-120.1200.1510.3430.2570.22320021-06-16
ISR0.8550.0500.0420.12832.5180.7270.0020-08-2020-12-180.0350.1000.3300.1420.1276221-09-07
ITA0.8730.2890.2520.00220.9820.8710.1220-05-2720-12-260.1140.2330.3290.3150.30618921-11-04
MEX0.7120.0380.0270.00440.0520.7070.2720-07-1320-12-230.1230.2340.5860.5590.53512422-11-23
NLD0.9290.3970.3690.00222.0390.9270.0720-06-1121-01-150.0720.1490.2010.1860.1798921-11-18
RUS0.9330.3370.3140.00080.4230.9330.0720-07-2220-12-140.0030.0490.1350.1180.1077421-10-30
SWE0.9220.6520.6010.00100.5400.9210.0820-10-1120-12-260.1260.1670.2600.2400.22916221-06-04
USA0.8680.2180.1890.00810.9480.8600.1220-09-0320-12-190.0940.1670.3260.2630.23815621-07-28
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Schlickeiser, R.; Kröger, M. Analytical Modeling of the Temporal Evolution of Epidemics Outbreaks Accounting for Vaccinations. Physics 2021, 3, 386-426. https://0-doi-org.brum.beds.ac.uk/10.3390/physics3020028

AMA Style

Schlickeiser R, Kröger M. Analytical Modeling of the Temporal Evolution of Epidemics Outbreaks Accounting for Vaccinations. Physics. 2021; 3(2):386-426. https://0-doi-org.brum.beds.ac.uk/10.3390/physics3020028

Chicago/Turabian Style

Schlickeiser, Reinhard, and Martin Kröger. 2021. "Analytical Modeling of the Temporal Evolution of Epidemics Outbreaks Accounting for Vaccinations" Physics 3, no. 2: 386-426. https://0-doi-org.brum.beds.ac.uk/10.3390/physics3020028

Article Metrics

Back to TopTop