Next Article in Journal
Large-Scale Simulation of Full Three-Dimensional Flow and Combustion of an Aero-Turbofan Engine on Sunway TaihuLight Supercomputer
Next Article in Special Issue
Unpredictable and Poisson Stable Oscillations of Inertial Neural Networks with Generalized Piecewise Constant Argument
Previous Article in Journal
Fair Numerical Algorithm of Coset Cardinality Spectrum for Distributed Arithmetic Coding
Previous Article in Special Issue
Influence of Information Blocking on the Spread of Virus in Multilayer Networks
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Epidemic Waves and Exact Solutions of a Sequence of Nonlinear Differential Equations Connected to the SIR Model of Epidemics

by
Nikolay K. Vitanov
1,2,* and
Kaloyan N. Vitanov
1
1
Institute of Mechanics, Bulgarian Academy of Sciences, Acad. G. Bonchev Str., Bl. 4, 1113 Sofia, Bulgaria
2
Climate, Atmosphere and Water Research Institute, Bulgarian Academy of Sciences, Blvd. Tzarigradsko Chaussee 66, 1784 Sofia, Bulgaria
*
Author to whom correspondence should be addressed.
Entropy 2023, 25(3), 438; https://doi.org/10.3390/e25030438
Submission received: 21 January 2023 / Revised: 21 February 2023 / Accepted: 27 February 2023 / Published: 1 March 2023

Abstract

:
The SIR model of epidemic spreading can be reduced to a nonlinear differential equation with an exponential nonlinearity. This differential equation can be approximated by a sequence of nonlinear differential equations with polynomial nonlinearities. The equations from the obtained sequence are treated by the Simple Equations Method (SEsM). This allows us to obtain exact solutions to some of these equations. We discuss several of these solutions. Some (but not all) of the obtained exact solutions can be used for the description of the evolution of epidemic waves. We discuss this connection. In addition, we use two of the obtained solutions to study the evolution of two of the COVID-19 epidemic waves in Bulgaria by a comparison of the solutions with the available data for the infected individuals.

1. Introduction

There are many complex systems of various scales around us. Examples range from atomic chains and lattices to systems of animals, humans and groups of humans, for example, research groups and social networks, economic systems, etc. [1,2,3,4,5,6,7]. These complex systems are usually nonlinear [8,9,10], and this complicates the study. One has to use the methodology of nonlinear time series analysis, and many of the corresponding models are based on differential or difference equations [11,12,13,14]. Then, one has to obtain solutions to these nonlinear equations. Usually, numerical methods are used to obtain the needed solutions. However, it is very useful if one can obtain exact analytical solutions to the model equations. In this case, one can easily study the relationships among the parameters of the complex system of interest. In addition, the exact solutions can be used as test solutions for checking the correctness of the work of the computer programs that have to supply numerical solutions to the corresponding systems of equations. In this article, we consider a nonlinear model of the evolution of epidemic waves (the SIR model of epidemics) and discuss a methodology for obtaining analytical solutions connected to this model. The methodology is based on the reduction of the model to a sequence of nonlinear differential equations. Several exact solutions to these nonlinear differential equations are obtained. Some of the solutions can be used for descriptions of epidemic waves.
Because of its importance, there is a large amount of research on the methodology for obtaining exact solutions to nonlinear differential equations. In the beginning, one tried to remove the nonlinearity of the studied equation by an appropriate transformation, such as the Hopf–Cole transformation [15,16]. Another transformation connected the nonlinear Korteweg–de Vries equation to the linear Schrödinger equation, and this led to the method of inverse scattering transform [17]. Other suitable transformations can be supplied by the truncated Painleve expansions [18,19,20,21,22]. Such a study led Kudryashov [23] to the formulation of the Method of Simplest Equation (MSE). The method is based on the determination of singularity order n of the solved NPDE. Then, a particular solution to this equation is searched as power series of solutions to a simpler equation called the simplest equation. For further results on this methodology and its applications, see [24,25,26,27,28].
Below we will use a methodology called SEsM (Simple Equations Method) [29,30,31,32,33]. Some specific cases of this methodology can be seen in our articles written approximately 30 years ago [34,35]. More than a decade ago [36], we used the ODE of Bernoulli as the simplest equation in the first version of the method: Modified Method of Simplest Equation (MMSE). We applied the methodology to population dynamics and ecology [37]. The MMSE [38] uses the concept of balanced equations for the fixation of the simplest equation. After this, the searched solution was constructed as a truncated power series of the solution to the simplest equation. This methodology leads to equivalent results with respect to the Method of Simplest Equation of Kudryashov.
We used MMSE actively till 2018 [39,40]. Our efforts have been constantly directed to extend the capacity of the methodology. The current version of the methodology (SEsM) can use more than one simple equation. SEsM based on two simple equations was used in [41]. The first discussion of SEsM was in [29]. Further discussion on SEsM can be seen in [33]. Applications of specific cases of SEsM can be seen in [42,43].
Below we apply SEsM to a system of nonlinear differential equations connected to an epidemic model: the SIR model of the spreading of epidemics in a population. There exist many models for the spread of epidemics. One of the most basic of these models is the SIR model for describing the temporal dynamics of an infectious disease in a population. The model compartmentalizes people into one of three categories: (i) susceptible to the disease; (ii) those who are currently infectious, and (iii) those who have recovered (with immunity). The SIR model is a set of equations that describes the number of people in each compartment at every point in time. A large amount of literature is devoted to this topic (for several examples, see [44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62]). Epidemic models can also be applied for the description of other processes, such as the spread of ideas (for overviews, see [3,63]). We also note the use of epidemic models for the study of COVID-19 spreading [64,65,66,67,68,69,70,71,72,73,74,75,76,77,78], as well as the numerical methods for obtaining solutions to such models [79,80].
We note the implicit solutions to the SIR model obtained by means of the Lambert function [81] and in [82]. Note also that the SIR model does not pass the Painleve test [81]. The solutions obtained below are explicit ones, and they are constructed by the use of elementary functions.
The text is organized as follows. In Section 2, we describe the methodology of SEsM. In Section 3, we obtain several exact solutions to the chain of nonlinear differential equations connected to the SIR model of epidemic spreading. Some of these solutions are appropriate to model epidemic waves. We discuss them in Section 4. In Section 5, we study the influence of the parameters of the SIR model on the shape of the epidemic waves described by the obtained solutions. In the same section, we use the obtained solutions to study two of the COVID-19 epidemic waves in Bulgaria. Several concluding remarks are summarized in Section 6.

2. Simple Equations Method (SEsM)

In general, the SEsM is constructed for obtaining exact and approximate solutions to systems of nonlinear differential equations. We will introduce a notation for the classification of the different cases of SEsM. In general, we have to solve n nonlinear differential equations. The main idea of SEsM is to obtain a solution on the basis of known solutions to m simpler differential equations. We denote this version of SEsM by SEsM(n,m). The most used version of SEsM up to now is this one for n = 1 . In this case, one has to solve one complicated nonlinear differential equation by means of the known solutions to m simpler differential equations. This will be denoted as SEsM(1,m). The most used specific version of SEsM(1,m) is the version SEsM(1,1). In this case, we solve one complicated nonlinear differential equation by means of the known solutions to one simpler differential equation. SEsM(1,1) is called also the Modified Method of Simplest Equation, and there are many applications of this specific case of SEsM [83,84,85].
The main idea of SEsM is as follows. We have to solve a system of nonlinear differential equations:
D i [ u i 1 ( x , , t ) , , u i n ( x , , t ) ] = 0 , i = 1 , 2 , , n .
D i [ u i 1 ( x , , t ) , , u i n ( x , , t ) , ] depends on the functions u i 1 ( x , , t ) , , u i n ( x , , t ) and some of their derivatives. The functions u i k can depend on several spatial coordinates. We have to transform (1) to
i = 1 n a i j ( ) E i j = 0 , j = 1 , 2 , , p i .
This happens by the representation of u i k as composite functions of known analytical solutions to simpler equations. Above, E i j are functions of the independent spatial variables and of time. The quantities a i j are very important for SEsM. a i j are relationships among the parameters of the Equation (1), parameters of the solutions and the parameters of the solutions to the simpler equations. p i is a parameter that is characteristic of the i-th equation from (1). It is important that the relationships a i j contain only parameters, whereas the spatial coordinates and the time are concentrated on the functions E i j . If we manage to reduce the Equation (1) to the form (2), then we can set
a i j ( ) = 0 .
Thus, we obtain a system of nonlinear algebraic equations. This system contains relationships among the parameters of (1), parameters of the used simpler equations and the parameters of the solution constructed by the simpler equations. Each nontrivial solution to (3) leads to a solution to the system (1).
SEsM has four steps. The first step of SEsM is connected to the possibility of transformations of the nonlinearities of the Equation (1). The solved complicated differential equations contain nonlinear combinations of the unknown functions and their derivatives. The experience shows that polynomial nonlinearity is the most treatable kind of nonlinearity in a nonlinear differential equation. Thus, if the nonlinearities in (1) are polynomial ones, then one does not need a transformation to convert these nonlinearities. However, if the nonlinearities are not polynomial ones, then one can use the transformation order to convert the nonlinearities to more treatable kinds of nonlinearities (and eventually to polynomial nonlinearities). Thus, one applies the transformations
u i k ( x , . . . , t ) = T i k [ F i k l ( x , , t ) , G i k l ( x , , t ) , ] .
In (4), T i k ( F i k l , G i k l , ) is a composite function of other functions F i k l , G i k l , . These other functions can be functions of several spatial variables and time. The transformations have two goals. First of all, and if it is possible, the transformations may remove the nonlinearity of the solved Equation (1). One example of this is the Hopf–Cole transformation. This transformation reduces the nonlinear Burgers equation to the linear heat equation [15,16]. However, removing the nonlinearity of an equation by means of a transformation is rarely achieved. Usually, the transformations T i k intend to transform the nonlinearity of the solved equations to a more treatable kind of nonlinearity, such as polynomial nonlinearity. For the specific case of SEsM(1,1), examples of the transformation T are: u ( x , t ) = 4 tanh 1 [ F ( x , t ) ] for the Poisson–Boltzmann equation, and u ( x , t ) = 4 tan 1 [ F ( x , t ) ] for the Sine–Gordon equation [34,35]. The Painleve expansion is another appropriate transformation. Finally, u ( x , . . . , t ) = F ( x , . . . , t ) (no transformation) is a possibility for certain classes of nonlinear differential equations. The application of (4) to (1) may lead to treatable nonlinear differential equations for F i k l , G i , k l , . The transformations T i k may remain unfixed at the first step of SEsM. Then, the functions T i k remain unknown, and we have to determine them at some of the following three steps of the methodology.
Step 2 of SEsM is connected to the construction of the solutions to the transformed Equation (1). The transformed equations contain the unknown functions and their derivatives. The basic idea of SEsM is to search for solutions to these equations in the form of composite functions containing solutions to simpler differential equations. The implementation of this idea leads to the necessity of working with derivatives of composite functions. These derivatives have to be expressed by the solutions to simpler equations and derivatives of these solutions. Because of this, one has to use the Faa di Bruno formula for the derivatives of the composite functions. In such a way, the solved equations are converted to relationships that contain functions that are solutions to more simple equations. At this point of the application of SEsM, the composite functions and the solutions to the more simple equation are still not fixed. We can fix the forms of the composite function in Step 2 of SEsM. However, in general, it is not necessary to do this at this step. One example of a fixation of the composite function for the case of SEsM(1,1) is
F = α + i 1 = 1 N β i 1 f i 1 + i 1 = 1 N i 2 = 1 N γ i 1 , i 2 f i 1 f i 2 + i 1 = 1 N i N = 1 N σ i 1 , , i N f i 1 f i N .
where α , β i 1 , γ i 1 , i 2 , σ i 1 , , i N are parameters, and f i k are functions that are solutions to more simple equations. The relationship used by Hirota [86] is a specific case of (5).
Step 3 of SEsM is usually the most important step in the application of the methodology. Here, we determine the form of the more simple equations of which the solutions will be used for solutions to the system of solved nonlinear differential equations. The simple equations, as well as the composite functions constructed by their solutions, are chosen in such a way that the solved differential equations are transformed into a system of relationships (1). One has to ensure that the relationships for a i j contain more than one term. Because of this, additional relationships among the parameters participating in the relationships for a i j may occur. These additional relationships are called balance equations.
At step 4 of SEsM, we use (3). This leads to a system of nonlinear algebraic relationships among the parameters of the Equation (1), the parameters of the composite functions, and the parameters of the solutions to the simpler equations. Any nontrivial solution to (3) leads to a solution to the system of the solved nonlinear equations.
For a specific case of applications of SEsM, see [29,30,31,32,33,36,37,38,39,40,41].

3. SEsM and Exact Analytical Solutions for a Chain of Equations Connected to the SIR Model of Epidemics

Below we use SEsM to obtain exact solutions to a chain of nonlinear differential equations connected to a specific nonlinear differential equation. This equation will be obtained on the basis of the SIR model from the epidemiology. Some of the obtained solutions can be used to the description of epidemic waves caused by different diseases (COVID-19 inclusive). Other obtained solutions will not be appropriate for this purpose.
The basic nonlinear differential equation is obtained by the use of the classic idea of Kermack and McKendrick [87] for the transformation of the SIR model with constant coefficients to a single nonlinear differential equation. We consider an epidemic in a population. The population is divided into three groups: susceptible individuals—S; infected individuals—I; recovered individuals—R. The model equations for the time change in the numbers of individuals from the above three groups are:
d S d t = τ N S I d I d t = τ N S I ρ I d R d t = ρ I .
In (6), τ is the transmission rate, and ρ is the recovery rate. These rates are assumed to be constants. From (6), we have the relationship
N = S + I + R .
N is the total population, which is assumed to be constant. The model (6) can be reduced to a single equation for R, as follows. From the last equation of (6), we have
I = 1 ρ d R d t .
The substitution of (8) in the first equation of (6) leads to the relationship
S = S ( 0 ) exp τ ρ N [ R R ( 0 ) ] .
Here, S ( 0 ) and R ( 0 ) are the numbers of susceptible individuals and those recovered at time t = 0 . The substitution of (7) and (9) in the last equation of (6) leads to the differential equation for R
d R d t = ρ N R S ( 0 ) exp τ ρ N ( R R ( 0 ) )
Below, we assume R ( 0 ) = 0 (no recovered individuals at t = 0 ). Let us consider the ratio τ R ρ N . We assume that τ R ρ N < < 1 . This can be realized, for example, when τ > ρ and R < < N . This means that we have an epidemic wave that affects a small amount of the population, and the number of recovered people remains small with respect to the number of the entire population. In this case, exp τ ρ N R can be represented as a Taylor series
exp τ ρ N R = j = 0 M τ ρ N R j .
M has infinite value in the full Taylor series, but we can truncate it at M = 2 , M = 3 ,..., if τ ρ N R is small enough. From (10), we obtain
d R d t = ρ N R S ( 0 ) j = 0 M τ ρ N R j , M = 2 , 3 ,
We set
α 0 = ρ [ N S ( 0 ) ] ; α 1 = τ S ( 0 ) N ρ ; α j = ( 1 ) j τ j S ( 0 ) ρ j 1 N j , j = 2 , 3 ,
Then (12) becomes
d R d t = j = 0 M α j R j
The chains of Equation (12) and (14) are connected to the orders of approximation of (10), which is the equation for the time evolution of the recovering individuals for an epidemic wave within the scope of the SIR model.
In (14), the independent variable is the time t. In principle, the independent variable can also be a spatial coordinate or a combination of spatial variables and time. In order to discuss this (more general) case below, we use an independent variable denoted as x. This variable can be a spatial variable, time, or some combination of spatial variables and time. Thus, we apply SEsM to the equation below
d R d x = j = 0 M α j R j ,
and we use the differential equations of Bernoulli and Riccati as simple equations. Our plan is as follows. First, we describe exact analytical solutions to the equations of the chain (14). Then, we adapt these solutions for the specific case of epidemics described by the chain of Equation (12).
First of all, we use the equation of Bernoulli
d y d x = p y + q y m , m = 2 , 3 , ,
as a simple equation within the scope of the SEsM methodology. By means of the transformation, y = u 1 / ( 1 m ) (16) is reduced to a linear differential equation. This leads to the solution to the equation of Bernoulli as follows
y ( x ) = p q + C p exp [ ( m 1 ) p x ] 1 m 1 .
In (17), C is a constant of integration.
We skip Step 1 of SEsM. No transformation is needed because the kind of nonlinearity of (15) is a polynomial one. In Step 2 of SEsM, we prescribe the composite function R ( y ) to be of the kind
R ( y ) = l = 0 L β l y l ,
where y ( x ) is the solution to the equation of Bernoulli and R ( y ) is the solution to (15). At Step 3 of SEsM, we have to obtain the balance equation. The presence of (16) and (18) fixes the balance equation of (15) to
m = 1 + L ( M 1 ) .
Thus, a specific solution to (15) has the form
R ( x ) = l = 0 L β l p q + C p exp { [ L ( M 1 ) ] p x } l L ( M 1 ) .
The parameters β l , p, q and C are fixed by the solution to the system of nonlinear algebraic equations at Step 4 of SEsM.
There is a specific case where we can obtain the general solution of (15). This case is M = 2 . Here, (15) becomes
d R d x = α 2 R 2 + α 1 R + α 0 .
(21) is an equation of the Riccati kind. For this equation, we know the specific solution
R ( x ) = α 1 2 α 2 θ 2 α 2 tanh θ ( x + C ) 2 ,
where θ 2 = α 1 2 4 α 0 α 2 > 0 , and C is a constant of integration. On the basis of the specific solution (22) of (21), we can write the general solution to (21) as R = α 1 2 α 2 θ 2 α 2 tanh θ ( x + C ) 2 + D v where D is a constant, and v ( t ) is the solution to the linear differential equation
d v d x θ tanh θ ( x + C ) 2 v = α 2 D
The solution to (23) is
v = cosh 2 θ ( x + C 2 E 2 α 2 D θ tanh θ ( x + C 2 ,
where E is a constant of integration. Then, the general solution to Equation (21) is
R ( x ) = α 1 2 α 2 θ 2 α 2 tanh θ ( x + C ) 2 + D cosh 2 θ ( x + C ) 2 E 2 α 2 D θ tanh θ ( x + C ) 2 .
Let us now obtain the form of several solutions to the kind (20). For M = 2 , we have the general solution (25) to the corresponding Equation (21). Thus, we start from M = 3 . The equation we have to solve is
d R d x = α 3 R 3 + α 2 R 2 + α 1 R + α 0 .
The solution is of the kind (18), and from (19), we have the balanced equation m = 1 + 2 L . This fixes the form of the simple equation of Bernoulli for this case: d y d x = p y + q y 1 + 2 L . We start from the simplest case L = 2 . The equation of Bernoulli becomes d y d x = p y + q y 5 , and the solution to (26) has the form R = β 2 y 2 + β 1 y + β 0 . The substitution of the last relationships in (26) leads to the following system of nonlinear algebraic relationships
2 q α 3 β 2 2 = 0 β 1 ( q 3 α 3 β 2 2 ) = 0 α 3 [ β 0 β 2 2 + 2 β 1 2 β 2 + β 2 ( 2 β 0 β 2 + β 1 2 ) ] α 2 β 2 2 = 0 β 1 [ 2 α 2 β 2 α 3 ( 6 β 0 β 2 + β 1 2 ) ] = 0 α 3 [ β 0 ( 2 β 0 β 2 + β 1 2 ) + 2 β 0 β 1 2 + β 0 2 β 2 ] α 2 ( 2 β 0 β 2 + β 1 2 ) α 1 β 2 = 0 β 1 ( p 3 α 3 β 0 2 2 α 2 β 0 α 1 ) = 0 α 3 β 0 3 + α 0 + α 1 β 0 + α 2 β 0 2 = 0
The solution to (27) is
q = α 3 β 2 2 2 ; β 1 = 0 ; β 0 = α 2 3 α 3 ; α 1 = α 2 2 3 α 3 ; α 0 = α 2 3 27 α 3 2 .
Thus, the equation
d R d x = α 3 R 3 + α 2 R 2 + α 2 2 3 α 3 R + α 2 3 27 α 3 2 ,
has the specific exact analytical solution
R ( x ) = l = 0 2 β l p q + C p exp { 4 p x } l 4 = α 2 3 α 3 + β 2 p α 3 β 2 2 2 + C p exp { 4 p x } 1 2
Next, we consider the case M = 3 , L = 3 . The equation of Bernoulli becomes d y d x = p y + q y 7 , and the solution to (26) has the form R = β 3 y 3 + β 2 y 2 + β 1 y + β 0 . The substitution of the last relationships in (26) leads to the following system of nonlinear algebraic relationships
3 q α 3 β 3 2 = 0 β 2 ( 2 q 3 α 3 β 3 2 ) = 0 β 1 q α 3 β 3 ( β 1 β 3 + 2 β 2 2 + 2 β 1 β 3 + β 2 2 ) = 0 α 2 β 3 2 α 3 β 3 ( 2 β 1 β 2 + 3 β 0 β 3 + 2 β 1 β 2 ) α 3 β 2 ( 2 β 1 β 3 + β 2 2 ) = 0 α 3 [ 2 β 0 β 2 β 3 + β 1 ( 2 β 1 β 3 + β 2 2 ) + β 2 ( 2 β 0 β 3 + 2 β 1 β 2 ) + β 3 ( 2 β 0 β 2 + β 1 2 ) ] 2 α 2 β 2 β 3 = 0 α 3 [ β 0 ( 2 β 1 β 3 + β 2 2 ) + β 1 ( 2 β 0 β 3 + 2 β 1 β 2 ) + β 2 ( 2 β 0 β 2 + β 1 2 ) + 2 β 0 β 1 β 2 ] α 2 ( 2 β 1 β 3 + β 2 2 ) = 0 3 β 3 p α 3 [ β 1 ( 2 β 0 β 2 + β 1 2 ) + β 0 ( 2 β 0 β 3 + 2 β 1 β 2 ) + β 3 β 0 2 + 2 β 0 β 1 β 2 ] α 2 ( 2 β 0 β 3 + 2 β 1 β 2 ) α 1 β 3 = 0 α 3 [ β 0 ( 2 β 0 β 2 + β 1 2 ) + 2 β 0 β 1 2 + β 2 β 0 2 ] α 2 ( 2 β 0 β 2 + β 1 2 ) α 1 β 2 = 0 β 1 ( p 3 α 3 β 0 2 2 α 2 β 0 α 1 ) = 0 α 0 α 3 β 0 3 α 2 β 0 2 α 1 β 0 = 0
The solution to (31) is
q = α 2 3 α 3 ; p = α 2 2 + 3 α 1 α 3 9 α 3 ; β 0 = α 2 3 α 3 ; β 1 = β 2 = 0 ; α 0 = α 2 ( 2 α 2 2 + 9 α 1 α 3 ) 27 α 3 3
Thus, the equation
d R d x = α 3 R 3 + α 2 R 2 + α 1 R + α 2 ( 2 α 2 2 + 9 α 1 α 3 ) 27 α 3 3 ,
has the specific exact analytical solution
R ( x ) = l = 0 3 β l p q + C p exp { 6 p x } l 6 = α 2 3 α 3 + β 3 α 2 2 + 3 α 1 α 3 3 α 2 + C ( α 2 2 + 3 α 1 α 3 ) exp 6 α 2 2 + 3 α 1 α 3 9 α 3 x 1 2
Next, we consider the case M = 4 . We have to solve the equation
d R d x = α 4 R 4 + α 3 R 3 + α 2 R 2 + α 1 R + α 0 .
The solution is of the kind (18), and from (19), we have the balanced equation m = 1 + 3 L . This fixes the form of the simple equation of Bernoulli for this case: d y d x = p y + q y 1 + 3 L . We start from the simplest case L = 2 . The equation of Bernoulli becomes d y d x = p y + q y 7 , and the solution to (35) has the form R = β 2 y 2 + β 1 y + β 0 . The substitution of the last relationships in (35) leads to the following system of nonlinear algebraic relationships
β 2 ( 2 q α 4 β 2 3 ) = 0 β 1 ( q 4 α 4 β 2 3 ) = 0 β 2 2 [ α 3 β 2 α 4 ( 2 β 0 β 2 + β 1 2 ) + 4 β 1 2 ] = 0 β 1 [ 3 α 3 β 2 2 α 4 ( 4 β 0 β 2 2 + 4 β 2 ( 2 β 0 β 2 + β 1 2 ) ) ) ] = 0 α 4 [ 2 β 0 2 β 2 2 + 8 β 0 β 1 2 β 2 + ( 2 β 0 β 2 + β 1 2 ) 2 ] α 3 [ β 0 β 2 2 + 2 β 1 2 β 2 + β 2 ( 2 β 0 β 2 + β 1 2 ) ] α 2 β 2 2 = 0 β 1 { α 4 [ 4 β 0 2 β 2 + 4 β 0 ( 2 β 0 β 2 + β 1 2 ) ] α 3 [ 4 β 0 β 2 + 2 β 0 β 2 + β 1 2 ] 2 α 2 β 2 } = 0 α 4 β 0 2 ( 4 β 0 β 2 + 6 β 1 2 ) α 3 β 0 ( 2 β 0 β 2 + 3 β 1 2 + β 2 ) α 2 ( 2 β 0 β 2 + β 1 2 ) α 1 β 2 = 0 β 1 ( p 4 α 4 β 0 3 3 α 3 β 0 2 2 α 2 β 0 α 1 ) = 0 α 3 β 0 3 α 0 α 2 β 0 2 α 4 β 0 4 α 1 β 0 = 0
The solution to (36) is
q = α 4 β 2 3 4 ; β 0 = α 3 4 α 4 ; β 1 = 0 ; α 0 = α 3 4 256 α 4 3 ; α 1 = α 3 2 16 α 4 2 ; α 2 = α 3 2 8 α 4
Thus, the equation
d R d x = α 4 R 4 + α 3 R 3 + α 3 2 8 α 4 R 2 + α 3 2 16 α 4 2 R + α 3 4 256 α 4 3 ,
has the solution
R = β 2 y 2 + β 1 y + β 0 = α 3 4 α 4 + β 2 p α 4 β 2 3 4 + C p exp { 6 p x } 1 3
Next, we consider the case M = 5 . We have to solve the equation
d R d x = α 5 R 5 + α 4 R 4 + α 3 R 3 + α 2 R 2 + α 1 R + α 0 .
The solution is of the kind (18), and from (19), we have the balanced equation m = 1 + 4 L . This fixes the form of the simple equation of Bernoulli for this case: d y d x = p y + q y 1 + 4 L . We start from the simplest case L = 2 . The equation of Bernoulli becomes d y d x = p y + q y 9 , and the solution of (35) has the form R = β 2 y 2 + β 1 y + β 0 . The substitution of the last relationships in (39) leads to the following system of nonlinear algebraic relationships
β 2 ( 2 q α 5 β 2 4 ) = 0 β 1 ( q 5 α 5 β 2 4 ) = 0 α 5 β 2 [ 5 β 0 β 2 3 + 10 β 1 2 β 2 2 + 4 β 0 β 1 β 2 2 + 4 β 1 ( 2 β 0 β 2 + β 1 2 ) ] = 0 β 1 [ 4 α 4 β 2 3 α 5 ( 8 β 0 β 2 3 + 4 β 1 2 β 2 2 + 4 β 1 β 2 2 ) + β 2 ( 12 β 0 β 2 2 + 4 β 1 2 β 2 ) ] = 0 α 5 [ 2 β 0 β 2 2 ( 2 β 0 β 2 + β 1 2 + 2 β 1 2 ) + β 1 2 β 2 ( 12 β 0 β 2 + 4 β 1 2 ) + β 2 ( 2 β 0 2 β 2 2 + 16 β 0 β 1 2 β 2 + 2 ( 2 β 0 β 2 + β 1 2 ) 2 ] α 4 β 2 2 ( 4 β 0 β 2 + 6 β 1 2 ) α 3 β 2 3 = 0 α 5 [ 4 β 0 β 1 β 2 ( β 0 β 2 + ( 2 β 0 β 2 + β 1 2 ) 2 ) + β 1 ( 2 β 0 2 β 2 2 + 8 β 0 β 1 2 β 2 + ( 2 β 0 β 2 + β 1 2 ) 2 ) + 4 β 0 β 1 β 2 ( β 0 β 2 + ( 2 β 0 β 2 + β 1 2 ) ) ] 4 α 4 β 1 β 2 ( 3 β 0 β 2 + β 1 2 ) 3 α 3 β 1 β 2 2 = 0 α 5 [ β 0 ( 2 β 0 2 β 2 2 + 8 β 0 β 1 2 β 2 + ( 2 β 0 β 2 + β 1 2 ) 2 ) + 4 β 0 β 1 2 ( 3 β 0 β 2 + β 1 2 ) + 2 β 0 2 β 2 ( 2 β 0 β 2 + 3 β 1 2 ) ] 3 α 4 β 2 ( β 0 β 2 + β 1 2 ) 3 α 3 β 2 ( β 0 β 2 + β 1 2 ) α 2 β 2 2 = 0 α 5 [ 4 β 0 2 β 1 ( 3 β 0 β 2 + β 1 2 ) + β 1 β 0 2 ( 4 β 0 β 2 + 6 β 1 2 ) + 4 β 2 β 0 3 β 1 ] 4 α 4 β 0 β 1 ( 2 β 0 β 2 + β 1 2 ) α 3 β 1 ( 6 β 0 β 2 + β 1 2 ) 2 α 2 β 1 β 2 = 0 α 5 β 0 3 ( 5 β 0 β 2 + 10 β 1 2 ) α 4 β 0 2 ( 4 β 0 β 2 + 6 β 1 2 ) 2 α 3 β 0 ( β 0 β 2 + β 1 2 ) α 2 ( β 0 β 2 + β 1 2 ) α 1 β 2 = 0 β 1 p 5 α 5 β 0 4 β 1 4 α 4 β 0 3 β 1 3 α 3 β 0 2 β 1 2 α 2 β 0 β 1 α 1 β 1 = 0 α 0 α 2 β 0 2 α 4 β 0 4 α 3 β 0 3 α 5 β 0 5 α 1 β 0 = 0
The solution to (41) is
q = α 5 β 2 4 2 ; β 0 = α 4 5 α 5 ; β 1 = 0 ; α 0 = α 4 5 3125 α 5 4 ; α 1 = α 4 4 125 α 5 3 ; α 2 = 2 α 4 3 25 α 5 2 ; α 3 = 2 α 4 2 5 α 5
Thus, the equation
d R d x = α 5 R 5 + α 4 R 4 + 2 α 4 2 5 α 5 R 3 + 2 α 4 3 25 α 5 2 R 2 + α 4 4 125 α 5 3 R + α 4 5 3125 α 5 4 ,
has the specific solution
R = β 2 p α 5 β 2 4 2 + C p exp { 8 p x } 1 4 α 4 5 α 5
Obtaining the exact solutions to the chain of equations can be continued. Below we focus on the epidemic waves connected to the SIR model. The additional exact solutions to the chain of equations will be discussed elsewhere.

4. Discussion of the Obtained Exact Analytical Solutions to the Studied Chain of Equations from the Point of View of Modeling of Epidemic Waves

Above, we have obtained several exact solutions to equations that can be connected to the SIR model of epidemic waves. The solutions are of two classes: (i) solutions that can be used for the purposes of the SIR model and (ii) solutions that are not appropriate for the use of the purposes of the SIR model. We begin the discussion by considering the solutions that can be used for the purposes of the SIR model. There are two groups of these solutions. The first group contains solutions without relationships among parameters α i in the form of equalities. The second group contains solutions with relationships among parameters α i , which are in the form of equalities.
Above, we have obtained relationships for the quantity R ( x ) , where x was some coordinate, which, in particular, can be some spatial coordinate, time, or a combination of time and spatial coordinates. Below, we consider the specific case when the coordinate x is time t. Thus, we obtain solutions to the number of recovered people R ( t ) for the case of the SIR model. This allows us to calculate the time evolution of infected people I on the basis of (8): I = 1 ρ d R d t . Then, we can calculate the relative growth rate
σ ( t ) = 1 I d I d t .
It can be written as
σ ( t ) = ρ ( R n 1 ) .
In (46)
R n ( t ) = 1 + σ ( t ) ρ ,
is the time-varying effective reproduction number. We use R n for the effective reproduction number in order to distinguish it from R by which we denote the recovered people. (46) shows that there is a specific value R n = 1 . If R n < 1 , then σ ( t ) < 0 , and the relative growth rate is negative. This means that d I / d t is negative; in other words, the number of infected individuals decreases, and the epidemic shrinks. If R n > 1 , then σ ( t ) > 0 , and the relative growth rate is positive. This means that d I / d t is positive; in other words, the number of infected individuals increases, and the epidemic extends. Note that we use (47) instead of the exact relationship R n = τ S ρ N , as we work with an approximate solution of (10).
We stress here the following. Our basic assumption for reducing the SIR model to a chain of equations was τ R ρ N 2 < < 1 . In order to keep the assumption in order, we have to consider epidemic waves for which R < < N . This means that the epidemic wave has to affect a small amount of the entire population. If this is not the case, we have to solve the SIR model numerically. For the rest of this section, we assume that we are within the scope of the assumption τ R ρ N 2 < < 1 , which allows us to obtain analytical results.
We have analytical relationships for several epidemic waves. Thus, we can calculate their characteristics by means of the relationships mentioned above. We denote R n ( 0 ) = τ S ( 0 ) ρ N . For the calculation of S, we will use the approximate relationship that occurs from (9)
S ( t ) = S ( 0 ) 1 τ R ρ N
We start from the specific solution (22). Taking into account (13) and (46), we obtain
R ( t ) = S ( 0 ) [ R n ( 0 ) 1 ] 2 R n ( 0 ) 2 + δ S ( 0 ) [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 R n ( 0 ) 2 tanh { 1 2 δ ρ × [ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) ] 1 / 2 { t + 2 artanh R n ( 0 ) 1 δ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 } } ,
where δ = ± 1 . Equation (47) allows us to calculate the time evolution of infected individuals for this specific solution. From (8), we obtain
I ( t ) = 1 ρ d R d t = S ( 0 ) 4 R n ( 0 ) 2 [ R n ( 0 ) 1 ] 2 + 4 R n 2 ( 0 ) N S ( 0 ) S ( 0 ) sech 2 { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 t + 2 artanh R n ( 0 ) 1 δ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 }
From (49) and (48), we obtain
S ( t ) = S ( 0 ) { 1 τ ρ N S ( 0 ) [ R n ( 0 ) 1 ] 2 R n ( 0 ) 2 + δ S ( 0 ) [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 R n ( 0 ) 2 tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 t + 2 artanh R n ( 0 ) 1 δ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 } } .
This allows us to calculate σ ( t ) from (46)
σ ( t ) = 1 I d I d t = δ ρ tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 t + 2 artanh R n ( 0 ) 1 δ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 } } .
Then, from (47)
R n ( t ) = 1 δ tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 t + 2 artanh R n ( 0 ) 1 δ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 } } .
We remember that the above results are valid if τ R ρ N 2 < < 1 . In other words, we have satisfy
{ τ ρ N { S ( 0 ) [ R n ( 0 ) 1 ] 2 R n ( 0 ) 2 + δ S ( 0 ) [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 R n ( 0 ) 2 tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 t + 2 artanh R n ( 0 ) 1 δ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 } } } 2 < < 1
We can obtain the following estimation for this condition. The maximum value of the tanh function is 1. This, from (14), we obtain
τ 2 ρ 2 N 2 S ( 0 ) [ R n ( 0 ) 1 ] 2 R n ( 0 ) 2 + δ S ( 0 ) [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 R n ( 0 ) 2 2 < < 1 .
Next we calculate the quantities for the solution (25). The requirement R ( 0 ) = 0 leads to the determination of the constant C as follows:
C = 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ]
The solution (25) becomes
R ( t ) = S ( 0 ) [ R n ( 0 ) 1 ] 2 R n ( 0 ) 2 + δ S ( 0 ) [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 R n ( 0 ) 2 tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } + D / { cosh 2 { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } { E + 2 D R n ( 0 ) 2 δ S ( 0 ) ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } } } .
In (57) δ = ± 1 . Equation (57) allows us to calculate the other quantities connected to this solution as follows
I = 1 ρ d R d t = S ( 0 ) 4 R n ( 0 ) 2 [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) × sech 2 { 1 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } } δ D ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 / { cosh 3 { 1 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } { E + 2 D R n ( 0 ) 2 δ S ( 0 ) ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 tanh { 1 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } } } D 2 R n ( 0 ) 2 S ( 0 ) sech 2 { 1 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } } / { cosh 2 { 1 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } { E + 2 D R n ( 0 ) 2 δ S ( 0 ) ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 tanh 2 { 1 2 δ ρ ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t +
2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } } } .
Thus,
S ( t ) = S ( 0 ) { 1 τ ρ N { S ( 0 ) [ R n ( 0 ) 1 ] 2 R n ( 0 ) 2 + δ S ( 0 ) [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 R n ( 0 ) 2 tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } + D / { cosh 2 { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } { E + 2 D R n ( 0 ) 2 δ S ( 0 ) ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } } } } }
Moreover,
σ ( t ) = { ρ S ( 0 ) T 1 3 / 2 tanh ( T 2 ) [ 1 tanh 2 ( T 2 ) ] 4 R n ( 0 ) 2 + 3 D ρ T 1 sinh 2 ( T 2 ) 2 T 3 cosh 4 ( T 2 ) + 2 δ D 2 ρ R n ( 0 ) 2 T 1 1 / 2 sinh ( T 2 ) [ 1 tanh 2 ( T 2 ) ] S ( 0 ) T 3 2 cosh 2 ( T 2 ) ρ D T 1 2 T 3 cosh 2 ( T 2 ) + 2 ρ D 3 R n ( 0 ) 4 [ 1 tanh 2 ( T 2 ) ] S ( 0 ) 2 T 3 3 cosh 2 ( T 2 ) + δ ρ D 2 R n ( 0 ) 2 T 1 1 / 2 tanh ( T 2 ) [ 1 tanh 2 ( T 2 ) ] S ( 0 ) T 3 2 cosh 2 ( T 2 ) } / { S ( 0 ) T 1 [ 1 tanh 2 ( T 2 ) ] 4 R n ( 0 ) 2 δ D T 1 1 / 2 sinh ( T 2 ) T 3 cosh 3 ( T 2 ) D 2 R n ( 0 ) 2 [ 1 tanh 2 ( T 2 ) ] S ( 0 ) T 3 2 cosh 2 ( T 2 ) }
where
T 1 = [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 [ N S ( 0 ) ] S ( 0 ) , T 2 = 1 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 [ N S ( 0 ) ] S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 [ N S ( 0 ) ] S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } , T 3 = E + 2 R n ( 0 ) 2 D tanh ( T 2 ) δ S ( 0 ) T 1 1 / 2 .
Finally,
R n ( t ) = 1 + { S ( 0 ) T 1 3 / 2 tanh ( T 2 ) [ 1 tanh 2 ( T 2 ) ] 4 R n ( 0 ) 2 + 3 D T 1 sinh 2 ( T 2 ) 2 T 3 cosh 4 ( T 2 ) + 2 δ D 2 R n ( 0 ) 2 T 1 1 / 2 sinh ( T 2 ) [ 1 tanh 2 ( T 2 ) ] S ( 0 ) T 3 2 cosh 2 ( T 2 ) D T 1 2 T 3 cosh 2 ( T 2 ) + 2 D 3 R n ( 0 ) 4 [ 1 tanh 2 ( T 2 ) ] S ( 0 ) 2 T 3 3 cosh 2 ( T 2 ) + δ D 2 R n ( 0 ) 2 T 1 1 / 2 tanh ( T 2 ) [ 1 tanh 2 ( T 2 ) ] S ( 0 ) T 3 2 cosh 2 ( T 2 ) } / { S ( 0 ) T 1 [ 1 tanh 2 ( T 2 ) ] 4 R n ( 0 ) 2 δ D T 1 1 / 2 sinh ( T 2 ) T 3 cosh 3 ( T 2 ) D 2 R n ( 0 ) 2 [ 1 tanh 2 ( T 2 ) ] S ( 0 ) T 3 2 cosh 2 ( T 2 ) }
The above results are valid if τ R ρ N 2 < < 1 . This means that
τ 2 ρ 2 N 2 { S ( 0 ) [ R n ( 0 ) 1 ] 2 R n ( 0 ) 2 + δ S ( 0 ) [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 R n ( 0 ) 2 tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } + D / { cosh 2 { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } { E + 2 D R n ( 0 ) 2 δ S ( 0 ) ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 tanh { 1 2 δ ρ × ( R n ( 0 ) 1 ) 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 { t + 2 δ ρ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 × artanh δ [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) 1 / 2 2 D R n ( 0 ) 2 S ( 0 ) E [ R n ( 0 ) 1 ] E [ R n ( 0 ) 1 ] 2 + 4 R n ( 0 ) 2 N S ( 0 ) S ( 0 ) + 2 D R n ( 0 ) 2 S ( 0 ) [ R n ( 0 ) 1 ] } } } } . } 2 < < 1
Next, we discuss the solution (34). For this solution, we have a relationship (28) for α 0 . From the point of view of the SIR model, α 0 > 0 , α 2 < 0 and α 3 > 0 . In order to ensure that α 0 > 0 , we see from (28) that it is necessary that 9 α 1 α 3 2 α 2 2 < 0 This leads to the requirement
R n ( 0 ) < 9 / 7 .
We note that 9 / 7 is close to 1.3 , which was a characteristic empirical value for the strong spreading of the virus in the case of the COVID-19 pandemic in recent years. Taking into account condition (64), we proceed with solution (34). The relationship (28) among the parameters α i fixes one of the parameters of the SIR model. Let us choose to fix S ( 0 ) . Then, we obtain
S ( 0 ) = N 7 ρ 2 N 2 + 27 τ 3 + δ 49 ρ 4 N 4 + 378 ρ 2 N 2 τ 3 + 729 τ 6 972 τ 2 ρ 3 N 2 1 / 2 54 τ 3 ,
where δ = ± 1 . From (64), it follows that
7 ρ 2 N 2 + 27 τ 3 + δ 49 ρ 4 N 4 + 378 ρ 2 N 2 τ 3 + 729 τ 6 972 τ 2 ρ 3 N 2 1 / 2 54 ρ τ 2 < 9 / 7 .
Below, we will write the solutions to δ = 1 as an example.
Another condition is R ( 0 ) = 0 , which fixes the value of the constant C to
C = { 9 τ [ 7 τ β 3 2 ρ 2 N 2 + 27 τ 4 β 3 2 τ β 3 2 49 ρ 4 N 4 + 378 ρ 2 N 2 τ 3 + 729 τ 6 972 τ 2 ρ 3 N 2 1 / 2 81 τ 3 ρ β 3 2 + 9 ρ 3 N 3 ] } / { ρ 2 N 2 [ 7 ρ 2 N 2 + 27 τ 3 49 ρ 4 N 4 + 378 ρ 2 N 2 τ 3 + 729 τ 6 972 τ 2 ρ 3 N 2 1 / 2 81 ρ τ 2 ] } .
Let us define
T 1 = 49 ρ 4 N 4 + 378 ρ 2 N 2 τ 3 + 729 τ 6 972 τ 2 ρ 3 N 2 ; T 2 = 7 ρ 2 N 2 + 27 τ 3 T 1 1 / 2 T 3 = T 2 2 2916 τ 2 ρ 2 N 2 + T 2 T 2 54 τ 2 ρ 18 ρ 2 N 2
Then, the solution (34) becomes
R ( t ) = ρ N 3 τ + β 3 { T 3 / [ [ T 2 18 τ ρ N + 9 τ T 3 ( 7 τ β 3 2 ρ 2 N 2 + 27 τ 4 β 3 2 τ β 3 2 T 1 1 / 2 81 τ 3 ρ β 3 2 + 9 ρ 3 N 3 ) exp 36 ρ 2 N 2 T 3 T 2 t ] / ρ 2 N 2 ( 7 ρ 2 N 2 + 27 τ 3 T 1 1 / 2 81 ρ τ 2 ) ] } 1 / 2
Then,
I = 162 β 3 T 3 3 T 5 τ exp 36 ρ 2 N 2 T 3 T 2 t / { ρ T 2 T 4 T 2 18 τ ρ N + 9 τ T 3 T 5 exp 36 ρ 2 N 2 T 3 T 2 t ρ 2 N 2 T 4 2 × T 3 T 2 18 τ ρ N + 9 τ T 3 T 5 exp 36 ρ 2 N 2 T 3 T 2 t ρ 2 N 2 T 4 1 / 2 } .
Here, T 1 , 2 , 3 are as above and
T 4 = T 2 81 ρ τ 2 ; T 5 = 27 τ 4 β 3 2 + 7 τ β 3 2 ρ 2 N 2 τ β 3 2 T 1 1 / 2 81 τ 3 ρ β 3 2 + 9 ρ 3 N 3 .
Furthermore,
σ ( t ) = ρ ( T 3 / T 6 ) 1 / 2 T 2 T 4 T 6 2 [ 26244 β 3 τ 2 T 3 6 T 5 2 T 7 2 ρ T 2 2 T 4 2 T 6 4 ( T 3 / T 6 ) 3 / 2 + 104976 β 3 τ 2 T 3 5 T 5 2 T 7 2 ρ T 2 2 T 4 2 T 6 3 ( T 3 / T 6 ) 1 / 2 5832 β 3 ρ τ N 2 T 3 4 T 5 T 7 T 2 2 T 4 T 6 2 ( T 3 / T 6 ) 1 / 2 ] / 162 β 3 τ T 3 3 T 5 T 7 .
Here T 1 , T 2 , T 3 , T 4 , T 5 are as above and
T 6 = T 2 18 τ ρ N + 9 τ T 3 T 5 T 7 ρ 2 N 2 T 4 , T 7 = exp 36 ρ 2 N 2 T 3 T 2 t .
Finally,
R n ( t ) = 1 + ( T 3 / T 6 ) 1 / 2 T 2 T 4 T 6 2 [ 26244 β 3 τ 2 T 3 6 T 5 2 T 7 2 ρ T 2 2 T 4 2 T 6 4 ( T 3 / T 6 ) 3 / 2 + 104976 β 3 τ 2 T 3 5 T 5 2 T 7 2 ρ T 2 2 T 4 2 T 6 3 ( T 3 / T 6 ) 1 / 2 5832 β 3 ρ τ N 2 T 3 4 T 5 T 7 T 2 2 T 4 T 6 2 ( T 3 / T 6 ) 1 / 2 ] / 162 β 3 τ T 3 3 T 5 T 7 .
The above results are valid if τ R ρ N 3 < < 1 . This means that
τ 3 ρ 3 N 3 { ρ N 3 τ + β 3 { T 3 / [ [ T 2 18 τ ρ N + 9 τ T 3 ( 7 τ β 3 2 ρ 2 N 2 + 27 τ 4 β 3 2 τ β 3 2 T 1 1 / 2 81 τ 3 ρ β 3 2 + 9 ρ 3 N 3 ) exp 36 ρ 2 N 2 T 3 T 2 t ] / ρ 2 N 2 ( 7 ρ 2 N 2 + 27 τ 3 T 1 1 / 2 81 ρ τ 2 ) ] } 1 / 2 } 3 < < 1
There are other solutions of this kind. They require specific values of one or several parameters connected to the epidemic wave. These specific values required decrease the probability of realization of the corresponding wave. Because of this, we will discuss this kind of specific solution elsewhere.
Next, we briefly discuss the exact solutions obtained, which are not appropriate for the use of the purposes of the SIR model. These solutions are (30), (39), and (44). The problems are connected to the values of the parameters α i corresponding to the SIR model. Let us consider the solution (30). From (28), we have the requirement α 0 = α 2 3 27 α 3 2 . However, from (13), it follows that α 2 < 0 and α 3 > 0 . Thus, we have α 0 = ρ [ N S ( 0 ) ] < 0 and then N < S ( 0 ) . The last relationship is false from the point of view of the SIR model. Thus, (30) is a valid solution, but it cannot be used for the purposes of the SIR model.
The next solution is (39). Here, we have the same problem with the parameter α 0 , which has to be positive. However, α 0 = α 3 4 / ( 256 α 4 3 ) and α 3 > 0 , α 4 < 0 for the specific case of the SIR model. Thus, we can not use (39) to model epidemic waves within the SIR model.
Next, we consider solution (44). Here, we have the same problem with α 0 , α 4 and α 5 as for the last two solutions for the specific case of the parameters of the SIR model. Therefore, we cannot use (44) to model epidemic waves within the SIR model.

5. Epidemic Waves Based on Some of the Obtained Solutions

Let us consider the influence of the parameters of the SIR model on the spread of the epidemic wave. The study will be made on the basis of some of the solutions obtained above.
Figure 1 shows the influence of the recovery rate ρ of the SIR model on the shape of the epidemic wave for the case of the relationship (50) obtained on the basis of solution (49). The decrease in the recovery rate leads to a larger peak of the wave (larger value of the maximum number of infected individuals for the studied wave). In addition, the peak of the wave occurs earlier. The increase in the value of the recovery rate ρ leads to a decrease in the maximum number of infected individuals. In addition, the peak of the epidemic wave is postponed, as can be seen from curves 3 and 4 in Figure 1a. The same kind of dependence on the maximum and the shape of the epidemic wave on the recovery rate ρ is observed for the relationship (58) for the epidemic wave obtained on the basis of solution (57) to the equation connected to the SIR model. Thus, the influence of the recovery rate on the epidemic wave is that the increased recovery rate leads to a faster decrease in the number of infected individuals, and this slows the rise of the epidemic wave and decreases its height.
Figure 2 shows the influence of the transmission rate τ on the shape of the epidemic wave. The increase in the transmission rate for the case of relationship (50) obtained by solution (49) leads to an increase in the value of the maximum number of infected individuals for the wave. In addition, the wave rises faster, as can be seen from curves 1 and 2 of Figure 2a. The effect of the decrease in the transmission rate on the shape of the wave described by the relationship (58) obtained by solution (57) is shown in Figure 2b. The decrease in τ , in this case, leads to a smaller maximum of the epidemic wave, and the wave occurs later. Thus, the increase in the transmission rate leads to a faster occurrence of the epidemic wave and an increase in the maximum number of infected individuals for the corresponding epidemic wave.
Figure 3 shows the influence of the initial number S ( 0 ) of susceptible individuals on the shape of the epidemic wave. We note that at t = 0 , we assume R ( 0 ) = 0 and then N = S ( 0 ) + I ( 0 ) . Then, the decrease in S ( 0 ) means that we have a larger value of I ( 0 ) . In other words, the decrease in the initial number of susceptible individuals means that the epidemic wave starts with a larger initial number of infected individuals. The influence of S ( 0 ) on the shape of the wave described by (50) obtained on the basis of solution (49) is shown in Figure 3a. The decrease in the initial number of susceptible people (the increase in the initial number of infected individuals) leads to a faster rise of the epidemic wave and a larger value of the maximum number of infected individuals. The result of the influence of S ( 0 ) on the epidemic wave described by relationship (58) (obtained on the basis of solution (57)) is the same, as can be seen in Figure 3b. Then, the larger value of suspected individuals at t = 0 (the smaller cluster of infected individuals at t = 0 ) leads to a later occurrence of the epidemic wave and a decrease in its height.
The above results of the influence of the parameters ρ , τ and S ( 0 ) on the shape of the epidemic wave hint at a strategy for fighting the epidemic. One needs to detect the epidemic when the cluster of infected individuals is still small. Then one has to try to decrease the transmission rate and increase the recovery rate. This can lead to a later occurrence of the epidemic wave and a decrease in the height of this wave.
The following figures show the influence of the parameters of the SIR model on the effective reproduction number R n connected to the epidemic wave. In principle, at the beginning of the wave, R n is larger than 1, and at the end of the wave, R n is smaller than 1. Figure 4 shows the influence of the recovery rate ρ on the effective reproduction number R n . Figure 4a shows the situation for the case of relationship (53) obtained on the basis of solution (49). We see that the decrease in the value of ρ leads to an increase in the initial value of R n , which is followed by a large decrease in the value of the effective reproduction number in the course of the value (see curves 1 and 2 of Figure 4a). The increase in the value of ρ results in a smaller initial value of R n and a smaller decrease in its value in the course of a wave. Large enough values of ρ lead to values of R n , which are close to 1 and correspond to a slowly rising epidemic wave.
Figure 4b shows the situation for the relationship (62) obtained on the basis of the solution (58). Quantitatively, the situation is the same as above. The increase in the value of the recovery rate leads to a decrease in the initial value of the effective reproduction number and to a smaller decrease in the value of this number in the course of the wave.
Figure 5 shows the influence of the transmission rate τ on the evolution of the effective reproduction number R n in the course of an epidemic wave. Figure 5a shows the situation for the case of relationship (53) obtained on the basis of solution (49). In this case, the decrease in the transmission rate leads to a decrease in the initial value of the effective reproduction number R n and to a smaller interval of decrease in the value of R n in the course of the epidemic wave. The same situation can be observed in Figure 5b for the case of relationship (62) obtained on the basis of the solution (58). We note that an appropriate value of the transmission rate combined with the corresponding values of the other parameters can make the value of R n closer to 2 and become even larger than this value.
Figure 6 shows the influence of the initial number of susceptible individuals S ( 0 ) on the value of the effective reproduction number R n . Figure 6a shows the situation for the case of relationship (53) obtained on the basis of solution (49). The decrease in the initial number of susceptible individuals (which corresponds to a larger number of infected individuals at t = 0 ) leads to a faster decrease in the value of the effective reproduction number R n . The same result can be seen in Figure 6b for relationship (62) obtained on the basis of solution (58).
Finally, we will use solutions (50) and (58) to approximate real data from the COVID-19 pandemic in Bulgaria. The data for the infected individuals for the first approximately 1000 days of the pandemic are shown in Figure 7. There have been several large COVID-19 epidemic waves in Bulgaria (the population of which is approximately 6.8 million people). In this article, we show how the above analytic results can be related to the second and third COVID-19 waves.
Figure 7 shows that there are periodic drops in the number of cases on Saturdays and Sundays and increases in the number of cases on Mondays. In order to remove this effect, which exists because of the presence of holidays, below we will work with the 7-day averages of the data I i * = 1 7 j = i 3 i + 3 I j and with the 14-day average of the data: I i * = 1 14 j = i 6 i + 7 I j .
Figure 8 shows the second COVID-19 wave in Bulgaria (dotted line shows the 7-day-average data) and its best fit with solutions (50) (Figure 8a) and (58) (Figure 8b). We observe that the fit with solution (58) is better, especially in the beginning and end regions of the wave. This better fit is also observed in the other figures below.
Figure 9 shows the fit of the 14-day averages of the data from solutions (50) and (58). Again, the fit by (58) is better. This can be expected as (58) is a more general solution in comparison to (50). The 14-day data are smoother than the 7-day-average data, and because of this, the fit of the 14-day-average data is better that the fit of the 7-day-average data.
Figure 10 and Figure 11 show the third large COVID-19 wave in Bulgaria and the corresponding fits of the 7-day-average data and 14-day-average data.
On the basis of the COVID-19 data and their fits, we can obtain the parameters of the models and compare these parameters for the two studied COVID-19 waves in Bulgaria. The comparison of the values of ρ (the recovery rate) obtained by the fits of the data for COVID-19 spreading in Bulgaria shows that ρ was larger for the second large wave in comparison to the third large wave. In addition, the transmission rate τ for the second large wave was larger in comparison to the transmission rate for the third large wave. Thus, the result is that the version of the COVID-19 virus that was responsible for the second large COVID- wave in Bulgaria spread faster than the version of the virus that was responsible for the third wave. In addition, the recovery time of the second large wave was faster in comparison to the recovery time of the third large wave.

6. Concluding Remarks

In this article, we apply the Simple Equations method (SEsM) to a chain of nonlinear differential equations connected to the SIR model of epidemic spreading. We obtain three classes of solutions from the point of their applicability for the purpose of epidemic modeling. The first class of solutions can be applied to model epidemics without imposing additional restrictions containing equalities among the parameters of the SIR model. Such solutions are (22) and (25). The second class of obtained solutions is solutions that require additional relationships containing equalities among the parameters of the SIR model. Such a solution is (34). The third class of the obtained solutions is solutions to the corresponding equation of the chain of equations but solutions that can not be used for the purpose of modeling epidemic spread. We note that the obtained solutions are appropriate for the description of a single epidemic wave that affects some populations and leads to an infection of a relatively small percentage of the individuals of this population. The obtained analytical solution allows us to study the influence of the parameters of the model (such as transmission rate, recovery rate and initial number of susceptible individuals) on the shape and evolution of the epidemic wave. The results are that larger recovery rates, smaller transmission rates and a larger number of potentially affected individuals (small number of infected individuals at the beginning of the wave) lead to a slower advancement of the wave and a decrease in its amplitude.
The obtained solution from the second class of solutions is also quite interesting. However, for its practical realization, it is required that specific relationships among the parameters of the epidemics are presented. There is more than one solution in this class of solutions to the chain of the studied equations. We intend to study these solutions elsewhere.
Finally, the third class of solutions demonstrates the capacity of SEsM. The methodology leads to numerous exact solutions to various equations, and this has been demonstrated many times already. The obtained solutions of this class can not be used to model epidemic waves as they lead to some relationships among the parameters of the SIR model that correspond to unrealistic assumptions about these parameters. Nevertheless, the obtained solutions are solutions to the corresponding equations from the studied chain of equations and can be used in other models where the corresponding relationships among the parameters lead to acceptable assumptions about the model parameters.

Author Contributions

Conceptualization: N.K.V.; methodology, N.K.V.; software, K.N.V.; validation, N.K.V.; formal analysis, N.K.V. and K.N.V.; resources, N.K.V.; data curation, K.N.V.; writing—original draft preparation, N.K.V. and K.N.V.; writing—review and editing, N.K.V.; visualization, K.N.V.; supervision, N.K.V.; project administration, N.K.V.; funding acquisition, N.K.V. All 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.

Data Availability Statement

Data for COVID-19 registered cases in Bulgaria can be found at https://coronavirus.bg/.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Latora, V.; Nicosia, V.; Russo, G. Complex Networks. Principles, Methods, and Applications; Cambridge University Press: Cambridge, UK, 2017; ISBN 978-1-107-10318-4. [Google Scholar]
  2. Chian, A.C.-L. Complex Systems Approach to Economic Dynamics; Springer: Berlin, Germany, 2007; ISBN 978-3-540-39752-6. [Google Scholar]
  3. Vitanov, N.K. Science Dynamics and Research Production. Indicators, Indexes, Statistical Laws and Mathematical Models; Springer: Cham, Switzerland, 2016; ISBN 978-3-319-41629-8. [Google Scholar]
  4. Treiber, M.; Kesting, A. Traffic Flow Dynamics: Data, Models, and Simulation; Springer: Berlin, Germany, 2013; ISBN 978-3-642-32460-4. [Google Scholar]
  5. Kutner, R.; Ausloos, M.; Grech, D.; Di Matteo, T.; Schinckus, C.; Stanley, H.E. Manifesto for a Post-Pandemic Modeling. Physica A 2019, 516, 240–253. [Google Scholar] [CrossRef] [Green Version]
  6. Simon, J.H. The Economic Consequences of Immigration; The University of Michigan Press: Ann Arbor, MI, USA, 1999; ISBN 978-0472086160. [Google Scholar]
  7. Dimitrova, Z.I. Flows of Substances in Networks and Network Channels: Selected Results and Applications. Entropy 2022, 24, 1485. [Google Scholar] [CrossRef]
  8. Drazin, P.G. Nonlinear Systems; Cambridge University Press: Cambridge, UK, 1992; ISBN 0-521-40489-4. [Google Scholar]
  9. Dimitrova, Z.I. Numerical Investigation of Nonlinear Waves Connected to Blood Flow in an Elastic Tube with Variable Radius. J. Theor. Appl. Mech. 2015, 45, 79–92. [Google Scholar] [CrossRef] [Green Version]
  10. Ganji, D.D.; Sabzehmeidani, Y.; Sedighiamiri, A. Nonlinear Systems in Heat Transfer; Elsevier: Amsterdam, The Netherlands, 2018; ISBN 978-0-12-812024-8. [Google Scholar]
  11. Kantz, H.; Schreiber, T. Nonlinear Time Series Analysis; Cambridge University Press: Cambridge, UK, 2004; ISBN 978-0511755798. [Google Scholar]
  12. Verhulst, F. Nonlinear Differential Equations and Dynamical Systems; Springer: Berlin, Germany, 2006; ISBN 978-3-540-60934-6. [Google Scholar]
  13. Mills, T. Applied Time Series Analysis; Academic Press: London, UK, 2019; ISBN 978-012-813117-6. [Google Scholar]
  14. Grossberg, S. Nonlinear Neural Networks: Principles, Mechanisms, and Architectures. Neural Netw. 1981, 1, 17–61. [Google Scholar] [CrossRef] [Green Version]
  15. Hopf, E. The Partial Differential Equation: ut + uux = ϵuxx. Commun. Pure Appl. Math. 1950, 3, 201–230. [Google Scholar] [CrossRef]
  16. Cole, J.D. On a Quasi-Linear Parabolic Equation Occurring in Aerodynamics. Q. Appl. Math. 1951, 9, 225–236. [Google Scholar] [CrossRef] [Green Version]
  17. Ablowitz, M.J.; Clarkson, P.A. Solitons, Nonlinear Evolution Equations and Inverse Scattering; Cambridge University Press: Cambridge, UK, 1991; ISBN 978-0511623998. [Google Scholar]
  18. Tabor, M. Chaos and Integrability in Dynamical Systems; Wiley: New York, NY, USA, 1989; ISBN 978-0471827283. [Google Scholar]
  19. Carrielo, F.; Tabor, M. Similarity Reductions from Extended Painleve Expansions for Nonintegrable Evolution Equations. Physica D 1991, 53, 59–70. [Google Scholar] [CrossRef]
  20. Carrielo, F.; Tabor, M. Painleve Expansions for Nonintegrable Evolution Equations. Physica D 1989, 39, 77–94. [Google Scholar] [CrossRef]
  21. Weiss, J.; Tabor, M.; Carnevalle, G. The Painleve Property for Partial Differential Equations. J. Math. Phys. 1983, 24, 522–526. [Google Scholar] [CrossRef]
  22. Kudryashov, N.A. On Types of Nonlinear Nonintegrable Equations with Exact Solutions. Phys. Lett. A 1991, 155, 269–275. [Google Scholar] [CrossRef]
  23. Kudryashov, N.A. Simplest Equation Method to Look for Exact Solutions of Nonlinear Differential Equations. Chaos Solitons Fractals 2005, 24, 1217–1231. [Google Scholar] [CrossRef] [Green Version]
  24. Kudryashov, N.A.; Loguinova, N.B. Extended Simplest Equation Method for Nonlinear Differential Equations. Appl. Math. Comput. 2008, 205, 361–365. [Google Scholar] [CrossRef]
  25. Kudryashov, N.A. Exact Solitary Waves of the Fisher Equation. Phys. Lett. A 2005, 342, 99–106. [Google Scholar] [CrossRef]
  26. Kudryashov, N.A. One Method for Finding Exact Solutions of Nonlinear Differential Equations. Commun. Nonlinear Sci. Numer. Simul. 2012, 17, 2248–2253. [Google Scholar] [CrossRef] [Green Version]
  27. Kudryashov, N.A. Highly Dispersive Optical Solitons of the Generalized Nonlinear Eighth-Order Schrödinger Equation. Optik 2020, 206, 164335. [Google Scholar] [CrossRef]
  28. Kudryashov, N.A. Solitary waves of the generalized Sasa-Satsuma equation with arbitrary refractive index. Optik 2021, 232, 166540. [Google Scholar] [CrossRef]
  29. Vitanov, N.K.; Dimitrova, Z.I.; Vitanov, K.N. Simple Equations Method (SEsM): Algorithm, Connection with Hirota Method, Inverse Scattering Transform Method, and Several Other Methods. Entropy 2021, 23, 10. [Google Scholar] [CrossRef]
  30. Vitanov, N.K. Simple Equations Method (SEsM): Review and New Results. AIP Conf. Ser. 2022, 2459, 020003. [Google Scholar] [CrossRef]
  31. Vitanov, N.K.; Dimitrova, Z.I. Simple Equations Method and Non-linear Differential Equations with Non-polynomial Non-linearity. Entropy 2021, 23, 1624. [Google Scholar] [CrossRef]
  32. Vitanov, N.K.; Dimitrova, Z.I.; Vitanov, K.N. On the Use of Composite Functions in the Simple Equations Method to Obtain Exact Solutions of Nonlinear Differential Equations. Computation 2021, 9, 104. [Google Scholar] [CrossRef]
  33. Vitanov, N.K. Simple Equations Method (SEsM): An Affective Algorithm for Obtaining Exact Solutions of Nonlinear Differential Equations. Entropy 2022, 24, 1653. [Google Scholar] [CrossRef] [PubMed]
  34. Martinov, N.; Vitanov, N. On Some Solutions of the Two-Dimensional Sine-Gordon Equation. J. Phys. A Math. Gen. 1992, 25, L419–L426. [Google Scholar] [CrossRef]
  35. Vitanov, N.K. Breather and Soliton Wave Families for the Sine-Gordon Equation. Proc. R. Soc. Lond. A 1998, 454, 2409–2423. [Google Scholar] [CrossRef]
  36. Vitanov, N.K.; Jordanov, I.P.; Dimitrova, Z.I. On Nonlinear Population Waves. Appl. Math. Comput. 2009, 215, 2950–2964. [Google Scholar] [CrossRef]
  37. Vitanov, N.K.; Dimitrova, Z.I. Application of The Method of Simplest Equation for Obtaining Exact Traveling-Wave Solutions for Two Classes of Model PDEs from Ecology and Population Dynamics. Commun. Nonlinear Sci. Numer. Simul. 2010, 15, 2836–2845. [Google Scholar] [CrossRef]
  38. Vitanov, N.K. Modified Method of Simplest Equation: Powerful Tool for Obtaining Exact and Approximate Traveling-Wave Solutions of Nonlinear PDEs. Commun. Nonlinear Sci. Numer. Simul. 2011, 16, 1176–1185. [Google Scholar] [CrossRef]
  39. Vitanov, N.K. On Modified Method of Simplest Equation for Obtaining Exact and Approximate Solutions of Nonlinear PDEs: The Role of the Simplest Equation. Commun. Nonlinear Sci. Numer. Simul. 2011, 16, 4215–4231. [Google Scholar] [CrossRef]
  40. Vitanov, N.K.; Dimitrova, Z.I.; Vitanov, K.N. Modified Method of Simplest Equation for Obtaining Exact Analytical Solutions of Nonlinear Partial Differential Equations: Further Development of the Methodology with Applications. Appl. Math. Comput. 2015, 269, 363–378. [Google Scholar] [CrossRef] [Green Version]
  41. Vitanov, N.K.; Dimitrova, Z.I. Modified Method of Simplest Equation Applied to the Nonlinear Schrödinger Equation. J. Theor. Appl. Mech., Sofia 2018, 48, 59–68. [Google Scholar] [CrossRef] [Green Version]
  42. Vitanov, N.K. Simple Equations Method (SEsM) and Its Connection with the Inverse Scattering Transform Method. AIP Conf. Proc. 2021, 2321, 030035. [Google Scholar] [CrossRef]
  43. Vitanov, N.K.; Dimitrova, Z.I. Simple Equations Method (SEsM) and Its Particular Cases: Hirota Method. AIP Conf. Proc. 2021, 2321, 030036. [Google Scholar] [CrossRef]
  44. Brauer, F.; Castillo-Chavez, C.; Feng, Z. Mathematcal Models in Epidemiology; Springer: New York, NY, USA, 2019; ISBN 978-1-4939-9828-9. [Google Scholar]
  45. Diekmann, O.; Heesterbeek, H.; Britton, T. Mathematical Tools for Understanding Infectious Disease Dynamics; Princeton University Press: Princeton, NJ, USA, 2012; ISBN 978-0-6911-5539-5. [Google Scholar]
  46. Martcheva, M. An Introduction to Mathematical Epidemiology; Springer: New York, NY, USA, 2015; ISBN 978-14899-7612-3. [Google Scholar]
  47. Li, M.I. An Introduction to Mathematical Modeling of Infectious Diseases; Springer: Cham, Switzerland, 2018; ISBN 978-3-319-72122-4. [Google Scholar]
  48. Brauer, F. Mathematical Epidemiology: Past, Present and Future. Infect. Dis. Model. 2017, 2, 113–127. [Google Scholar] [CrossRef] [PubMed]
  49. Britton, T. Stochastic Epidemic Models: A Survey. Math. Biosci. 2010, 225, 24–35. [Google Scholar] [CrossRef] [PubMed]
  50. Hethcote, H.W. A Thousand and One Epidemic Models. In Frontiers in Mathematical Biology; Levin, S.A., Ed.; Springer: Berlin, Germany, 1994; pp. 504–515. ISBN 978-3-642-50126-5. [Google Scholar]
  51. Keeling, M.J.; Eames, K.T. Networks and Epidemic Models. J. R. Soc. Interface 2005, 2, 295–307. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Capasso, V.; Serio, G. A Generalization of the Kermack- McKendrick Deterministic Epidemic Model. Math. Biosci. 1978, 42, 43–61. [Google Scholar] [CrossRef]
  53. Teng, P.S. A Comparison of Simulation Approaches to Epidemic Modeling. Annu. Rev. Phytopathol. 1985, 23, 351–379. [Google Scholar] [CrossRef]
  54. Hethcote, H.W. The Mathematics of Infectious Diseases. SIAM Rev. 2000, 42, 599–653. [Google Scholar] [CrossRef] [Green Version]
  55. Wang, W.; Tang, M.; Stanley, H.E.; Braunstein, L.A. Unification of Theoretical Approaches for Epidemic Spreading on Complex Networks. Rep. Prog. Phys. 2017, 80, 036603. [Google Scholar] [CrossRef]
  56. Cifuentes-Faura, J.; Faura-Martínez, U.; Lafuente-Lechuga, M. Mathematical Modeling and the Use of Network Models as Epidemiological Tools. Mathematics 2022, 10, 3347. [Google Scholar] [CrossRef]
  57. Rahimi, I.; Gandomi, A.H.; Asteris, P.G.; Chen, F. Analysis and Prediction of COVID-19 Using SIR, SEIQR, and Machine Learning Models: Australia, Italy, and UK Cases. Information 2021, 12, 109. [Google Scholar] [CrossRef]
  58. Cui, Q.; Qiu, Z.; Liu, W.; Hu, Z. Complex Dynamics of an SIR Epidemic Model with Nonlinear Saturate Incidence and Recovery Rate. Entropy 2017, 19, 305. [Google Scholar] [CrossRef]
  59. Trawicki, M.B. Deterministic Seirs Epidemic Model for Modeling Vital Dynamics, Vaccinations, and Temporary Immunity. Mathematics 2017, 5, 7. [Google Scholar] [CrossRef] [Green Version]
  60. Kozioł, K.; Stanislawski, R.; Bialic, G. Fractional-Order SIR Epidemic Model for Transmission Prediction of COVID-19 Disease. Appl. Sci. 2020, 10, 8316. [Google Scholar] [CrossRef]
  61. Godio, A.; Pace, F.; Vergnano, A. SEIR Modeling of the Italian Epidemic of SARS-CoV-2 Using Computational Swarm Intelligence. Int. J. Environ. Res. Public Health 2020, 17, 3535. [Google Scholar] [CrossRef]
  62. Frank, T.D. COVID-19 Epidemiology and Virus Dynamics; Springer: Cham, Switzerland, 2022; ISBN 978-3-030-97178-6. [Google Scholar]
  63. Vitanov, N.K.; Ausloos, M.R. Knowledge Epidemics and Population Dynamics Models for Describing Idea Diffusion. In Models of Science Dynamics; Scharnhorst, A., Boerner, K., Besselaar, P., Eds.; Springer: Berlin, Germany, 2010; pp. 69–125. ISBN 978-3-642-23068-4. [Google Scholar]
  64. Al-Shbeil, I.; Djenina, N.; Jaradat, A.; Al-Husban, A.; Ouannas, A.; Grassi, G. A New COVID-19 Pandemic Model including the Compartment of Vaccinated Individuals: Global Stability of the Disease-Free Fixed Point. Mathematics 2023, 11, 576. [Google Scholar] [CrossRef]
  65. Lee, S.J.; Lee, S.E.; Kim, J.-O.; Kim, G.B. Two-Way Contact Network Modeling for Identifying the Route of COVID-19 Community Transmission. Informatics 2021, 8, 22. [Google Scholar] [CrossRef]
  66. Harjule, P.; Poonia, R.C.; Agrawal, B.; Saudagar, A.K.J.; Altameem, A.; Alkhathami, M.; Khan, M.B.; Hasanat, M.H.A.; Malik, K.M. An Effective Strategy and Mathematical Model to Predict the Sustainable Evolution of the Impact of the Pandemic Lockdown. Healthcare 2022, 10, 759. [Google Scholar] [CrossRef]
  67. Etxeberria-Etxaniz, M.; Alonso-Quesada, S.; De la Sen, M. On an SEIR Epidemic Model with Vaccination of Newborns and Periodic Impulsive Vaccination with Eventual On-Line Adapted Vaccination Strategies to the Varying Levels of the Susceptible Subpopulation. Appl. Sci. 2020, 10, 8296. [Google Scholar] [CrossRef]
  68. Nkague Nkamba, L.; Manga, T.T. Modelling and Prediction of the Spread of COVID-19 in Cameroon and Assessing the Governmental Measures (March–September 2020). COVID 2021, 1, 622–644. [Google Scholar] [CrossRef]
  69. Almeshal, A.M.; Almazrouee, A.I.; Alenizi, M.R.; Alhajeri, S.N. Forecasting the Spread of COVID-19 in Kuwait Using Compartmental and Logistic Regression Models. Appl. Sci. 2020, 10, 3402. [Google Scholar] [CrossRef]
  70. Chen, J.; Yin, T. Transmission Mechanism of Post-COVID-19 Emergency Supply Chain Based on Complex Network: An Improved SIR Model. Sustainability 2023, 15, 3059. [Google Scholar] [CrossRef]
  71. Batool, H.; Li, W.; Sun, Z. Extinction and Ergodic Stationary Distribution of COVID-19 Epidemic Model with Vaccination Effects. Symmetry 2023, 15, 285. [Google Scholar] [CrossRef]
  72. Khorev, V.; Kazantsev, V.; Hramov, A. Effect of Infection Hubs in District-Based Network Epidemic Spread Model. Appl. Sci. 2023, 13, 1194. [Google Scholar] [CrossRef]
  73. Jitsinchayakul, S.; Humphries, U.W.; Khan, A. The SQEIRP Mathematical Model for the COVID-19 Epidemic in Thailand. Axioms 2023, 12, 75. [Google Scholar] [CrossRef]
  74. Ni, G.; Wang, Y.; Gong, L.; Ban, J.; Li, Z. Parameters Sensitivity Analysis of COVID-19 Based on the SCEIR Prediction Model. COVID 2022, 2, 1787–1805. [Google Scholar] [CrossRef]
  75. Wang, W.; Xia, Z. Study of COVID-19 Epidemic Control Capability and Emergency Management Strategy Based on Optimized SEIR Model. Mathematics 2023, 11, 323. [Google Scholar] [CrossRef]
  76. Leonov, A.; Nagornov, O.; Tyuflin, S. Modeling of Mechanisms of Wave Formation for COVID-19 Epidemic. Mathematics 2023, 11, 167. [Google Scholar] [CrossRef]
  77. Margenov, S.; Popivanov, N.; Ugrinova, I.; Hristov, T. Mathematical Modeling and Short-Term Forecasting of the COVID-19 Epidemic in Bulgaria: SEIRS Model with Vaccination. Mathematics 2022, 10, 2570. [Google Scholar] [CrossRef]
  78. Chang, Y.-C.; Liu, C.-T. A Stochastic Multi-Strain SIR Model with Two-Dose Vaccination Rate. Mathematics 2022, 10, 1804. [Google Scholar] [CrossRef]
  79. Noeiaghdam, S.; Micula, S. Dynamical Strategy to Control the Accuracy of the Nonlinear Bio-Mathematical Model of Malaria Infection. Mathematics 2021, 9, 1031. [Google Scholar] [CrossRef]
  80. Noeiaghdam, S.; Micula, S.; Nieto, J.J. A Novel Technique to Control the Accuracy of a Nonlinear Fractional Order Model of COVID-19: Application of the CESTAC Method and the CADNA Library. Mathematics 2021, 9, 1321. [Google Scholar] [CrossRef]
  81. Kudryashov, N.A.; Chmykhov, M.A.; Vigdorowitsch, M. Analytical Features of the SIR Model and their Applications to COVID-19. Appl. Math. Model. 2021, 90, 466–473. [Google Scholar] [CrossRef] [PubMed]
  82. Harko, T.; Lobo, F.S.N.; Mak, M.K. Exact Analytical Solutions of the Susceptible-Infected-Recovered (SIR) Epidemic Model and of the SIR Model with Equal Death and Birth Rates. Appl. Math. Comput. 2014, 236, 184–194. [Google Scholar] [CrossRef] [Green Version]
  83. Dimitrova, Z.I. Relation Between G’/G-expansion Method and the Modified Method of Simplest Equation. C. R. L’Acad. Bulg. Des Sci. 2012, 65, 1513–1520. [Google Scholar]
  84. Dimitrova, Z. On Traveling Waves in Lattices: The Case of Riccati Lattices. J. Theor. Appl. Mech. 2012, 42, 3–22. [Google Scholar] [CrossRef] [Green Version]
  85. Dimitrova, Z.I. Several Examples of Application of the Simple Equations Method (SEsM) for Obtaining Exact Solutions of Nonlinear PDEs. AIP Conf. Proc. 2022, 2459, 030005. [Google Scholar] [CrossRef]
  86. Hirota, R. The Direct Method in Soliton Theory; Cambridge University Press: Cambridge, UK, 2004; ISBN 0-521-83660-3. [Google Scholar]
  87. Kermack, W.O.; McKendrick, A.G. A Contribution to the Mathematical Theory of Epidemics. Proc. R. Soc. Lond. Ser. A 1927, 115, 700–721. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Influence of the recovery rate ρ on the number of infected people. Figure (a): the solution (50). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of the parameter ρ . Curve 2: ρ = 0.0057 , Curve 3: ρ = 0.007 , Curve 4: ρ = 0.008 . Figure (b): the solution (58). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are changes only in the value of the parameter ρ . Curve 2: ρ = 0.003 = Curve 3: ρ = 0.004 . Curve 4: ρ = 0.005 . Curve 5: ρ = 0.0065 .
Figure 1. Influence of the recovery rate ρ on the number of infected people. Figure (a): the solution (50). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of the parameter ρ . Curve 2: ρ = 0.0057 , Curve 3: ρ = 0.007 , Curve 4: ρ = 0.008 . Figure (b): the solution (58). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are changes only in the value of the parameter ρ . Curve 2: ρ = 0.003 = Curve 3: ρ = 0.004 . Curve 4: ρ = 0.005 . Curve 5: ρ = 0.0065 .
Entropy 25 00438 g001
Figure 2. Influence of the transmission rate τ on the number of infected people. Figure (a): the solution (50). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of the parameter τ . Curve 2: τ = 0.0095 , Curve 3: τ = 0.008 , Curve 4: τ = 0.007 . Figure (b): the solution (58). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are changes only in the value of the parameter τ . Curve 2: τ = 0.0115 . Curve 3: τ = 0.0105 . Curve 4: ρ = 0.0095 . Curve 5: ρ = 0.0088 .
Figure 2. Influence of the transmission rate τ on the number of infected people. Figure (a): the solution (50). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of the parameter τ . Curve 2: τ = 0.0095 , Curve 3: τ = 0.008 , Curve 4: τ = 0.007 . Figure (b): the solution (58). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are changes only in the value of the parameter τ . Curve 2: τ = 0.0115 . Curve 3: τ = 0.0105 . Curve 4: ρ = 0.0095 . Curve 5: ρ = 0.0088 .
Entropy 25 00438 g002
Figure 3. Influence of the parameter S ( 0 ) on the number of infected people. Figure (a): the solution (50). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are only changes in the value of parameter S ( 0 ) . Curve 2: S ( 0 ) = 999,999, Curve 3: S ( 0 ) = 999,000, Curve 4: S ( 0 ) = 998,000. Figure (b): the solution (58). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are changes only in the value of parameter S ( 0 ) . Curve 2: S ( 0 ) = 999,500. Curve 3: S ( 0 ) = 998,500. Curve 4: S ( 0 ) = 995,000. Curve 5: S ( 0 ) = 990,000.
Figure 3. Influence of the parameter S ( 0 ) on the number of infected people. Figure (a): the solution (50). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are only changes in the value of parameter S ( 0 ) . Curve 2: S ( 0 ) = 999,999, Curve 3: S ( 0 ) = 999,000, Curve 4: S ( 0 ) = 998,000. Figure (b): the solution (58). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are changes only in the value of parameter S ( 0 ) . Curve 2: S ( 0 ) = 999,500. Curve 3: S ( 0 ) = 998,500. Curve 4: S ( 0 ) = 995,000. Curve 5: S ( 0 ) = 990,000.
Entropy 25 00438 g003
Figure 4. Influence of the recovery rate ρ on the effective reproduction number R n from Equation (47). Figure (a): the relationship (53) obtained on the basis of the solution (49). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of the parameter ρ . Curve 2: ρ = 0.0057 , Curve 3: ρ = 0.007 , Curve 4: ρ = 0.008 . Figure (b): the relationship (62) obtained on the basis of the solution (58). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are only changes in the value of the parameter ρ . Curve 2: ρ = 0.0065 . Curve 3: ρ = 0.005 .
Figure 4. Influence of the recovery rate ρ on the effective reproduction number R n from Equation (47). Figure (a): the relationship (53) obtained on the basis of the solution (49). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of the parameter ρ . Curve 2: ρ = 0.0057 , Curve 3: ρ = 0.007 , Curve 4: ρ = 0.008 . Figure (b): the relationship (62) obtained on the basis of the solution (58). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are only changes in the value of the parameter ρ . Curve 2: ρ = 0.0065 . Curve 3: ρ = 0.005 .
Entropy 25 00438 g004
Figure 5. Influence of the transmission rate τ on the effective reproduction number R n Equation (47). Figure (a): the relationship (53). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of the parameter τ . Curve 2: τ = 0.0095 , Curve 3: τ = 0.008 , Curve 4: τ = 0.007 . Figure (b): the relationship (62). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other cures, there are changes only in the value of the parameter τ . Curve 2: τ = 0.0105 . Curve 3: τ = 0.0115 .
Figure 5. Influence of the transmission rate τ on the effective reproduction number R n Equation (47). Figure (a): the relationship (53). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of the parameter τ . Curve 2: τ = 0.0095 , Curve 3: τ = 0.008 , Curve 4: τ = 0.007 . Figure (b): the relationship (62). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other cures, there are changes only in the value of the parameter τ . Curve 2: τ = 0.0105 . Curve 3: τ = 0.0115 .
Entropy 25 00438 g005
Figure 6. Influence of the parameter S ( 0 ) on the effective reproduction number R n Equation (47). Figure (a): the relationship (53). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of parameter S ( 0 ) . Curve 2: S ( 0 ) = 999,999, Curve 3: S ( 0 ) = 999,000, Curve 4: S ( 0 ) = 998,000. Figure (b): the relationship (62). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are changes only in the value of the parameter S ( 0 ) . Curve 2: S ( 0 ) = 999,500. Curve 3: S ( 0 ) = 998,500. Curve 4: S ( 0 ) = 995,000. Curve 5: S ( 0 ) = 990,000.
Figure 6. Influence of the parameter S ( 0 ) on the effective reproduction number R n Equation (47). Figure (a): the relationship (53). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 . For the other curves, there are changes only in the value of parameter S ( 0 ) . Curve 2: S ( 0 ) = 999,999, Curve 3: S ( 0 ) = 999,000, Curve 4: S ( 0 ) = 998,000. Figure (b): the relationship (62). Curve 1: basic solution with parameters N = 10 6 , S ( 0 ) = 999,990, τ = 0.009 , ρ = 0.006 , D = 10 5 , E = 1 . For the other curves, there are changes only in the value of the parameter S ( 0 ) . Curve 2: S ( 0 ) = 999,500. Curve 3: S ( 0 ) = 998,500. Curve 4: S ( 0 ) = 995,000. Curve 5: S ( 0 ) = 990,000.
Entropy 25 00438 g006
Figure 7. COVID-19 epidemic waves in Bulgaria. X-axis: days since the beginning of the pandemic in Bulgaria (8-th of March 2020). Y-axis: registered number I of infected people per day. Wave 2 and wave 3 will be compared to the analytical results in this article.
Figure 7. COVID-19 epidemic waves in Bulgaria. X-axis: days since the beginning of the pandemic in Bulgaria (8-th of March 2020). Y-axis: registered number I of infected people per day. Wave 2 and wave 3 will be compared to the analytical results in this article.
Entropy 25 00438 g007
Figure 8. The second large COVID-19 wave in Bulgaria and the best fit of the 7-day-averaged data with the solutions (50) and (58). Dots: infected people (7-day average). Solid curves: Figure (a): solution (50). ρ = 0.0000982 , τ = 0.007878 . Figure (b): solution (58). ρ = 0.0000893 , τ = 0.007883 .
Figure 8. The second large COVID-19 wave in Bulgaria and the best fit of the 7-day-averaged data with the solutions (50) and (58). Dots: infected people (7-day average). Solid curves: Figure (a): solution (50). ρ = 0.0000982 , τ = 0.007878 . Figure (b): solution (58). ρ = 0.0000893 , τ = 0.007883 .
Entropy 25 00438 g008
Figure 9. The second large COVID-19 wave in Bulgaria and the best fit of the 14-day-average data from the solutions (50) and (58). Dots: infected people (14-day average). Solid curves: Figure (a): solution (50). ρ = 0.0000985 , τ = 0.007875 . Figure (b): solutions (58). ρ = 0.0000813 , τ = 0.007863 .
Figure 9. The second large COVID-19 wave in Bulgaria and the best fit of the 14-day-average data from the solutions (50) and (58). Dots: infected people (14-day average). Solid curves: Figure (a): solution (50). ρ = 0.0000985 , τ = 0.007875 . Figure (b): solutions (58). ρ = 0.0000813 , τ = 0.007863 .
Entropy 25 00438 g009
Figure 10. The third large COVID-19 wave in Bulgaria and the best fit of the 7-day-averaged data by the solutions (50) and (58). Dots: infected people (7-day average). Solid curves: Figure (a): solution (50). ρ = 0.0000598 , τ = 0.005285 . Figure (b): solutions (58). ρ = 0.0000599 , τ = 0.005296 .
Figure 10. The third large COVID-19 wave in Bulgaria and the best fit of the 7-day-averaged data by the solutions (50) and (58). Dots: infected people (7-day average). Solid curves: Figure (a): solution (50). ρ = 0.0000598 , τ = 0.005285 . Figure (b): solutions (58). ρ = 0.0000599 , τ = 0.005296 .
Entropy 25 00438 g010
Figure 11. The third large COVID-19 wave in Bulgaria and the best fit of the 14-day-averaged data from solutions (50) and (58). Dots: infected people 14-day average). Solid curves: Figure (a): solution (50). ρ = 0.0000676 , τ = 0.005326 . Figure (b): solutions (58). ρ = 0.0000699 , τ = 0.005299 .
Figure 11. The third large COVID-19 wave in Bulgaria and the best fit of the 14-day-averaged data from solutions (50) and (58). Dots: infected people 14-day average). Solid curves: Figure (a): solution (50). ρ = 0.0000676 , τ = 0.005326 . Figure (b): solutions (58). ρ = 0.0000699 , τ = 0.005299 .
Entropy 25 00438 g011
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Vitanov, N.K.; Vitanov, K.N. Epidemic Waves and Exact Solutions of a Sequence of Nonlinear Differential Equations Connected to the SIR Model of Epidemics. Entropy 2023, 25, 438. https://doi.org/10.3390/e25030438

AMA Style

Vitanov NK, Vitanov KN. Epidemic Waves and Exact Solutions of a Sequence of Nonlinear Differential Equations Connected to the SIR Model of Epidemics. Entropy. 2023; 25(3):438. https://doi.org/10.3390/e25030438

Chicago/Turabian Style

Vitanov, Nikolay K., and Kaloyan N. Vitanov. 2023. "Epidemic Waves and Exact Solutions of a Sequence of Nonlinear Differential Equations Connected to the SIR Model of Epidemics" Entropy 25, no. 3: 438. https://doi.org/10.3390/e25030438

APA Style

Vitanov, N. K., & Vitanov, K. N. (2023). Epidemic Waves and Exact Solutions of a Sequence of Nonlinear Differential Equations Connected to the SIR Model of Epidemics. Entropy, 25(3), 438. https://doi.org/10.3390/e25030438

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop