==== Front
PLoS One
PLoS One
plos
PLOS ONE
1932-6203
Public Library of Science San Francisco, CA USA

PONE-D-23-03892
10.1371/journal.pone.0287556
Research Article
Physical Sciences
Mathematics
Differential Equations
Medicine and Health Sciences
Pharmacology
Pharmacologic Analysis
Pharmacokinetic Analysis
Compartment Models
Research and Analysis Methods
Mathematical and Statistical Techniques
Mathematical Functions
Computer and Information Sciences
Systems Science
Nonlinear Systems
Physical Sciences
Mathematics
Systems Science
Nonlinear Systems
Research and Analysis Methods
Simulation and Modeling
Computer and Information Sciences
Systems Science
Nonlinear Dynamics
Physical Sciences
Mathematics
Systems Science
Nonlinear Dynamics
Physical Sciences
Mathematics
Numerical Analysis
Interpolation
Medicine and Health Sciences
Medical Conditions
Infectious Diseases
Infectious Disease Control
Direct numerical solutions of the SIR and SEIR models via the Dirichlet series approach
Direct numerical solutions of the SIR and SEIR models via the Dirichlet series approach
https://orcid.org/0000-0002-0463-8162
Prathom Kiattisak Conceptualization Data curation Formal analysis Investigation Methodology Software Validation Visualization Writing – original draft Writing – review & editing *
Jampeepan Asama Formal analysis Software
Division of Mathematics and Statistics, Walailak University, Nakhon Si Thammarat, Thailand
Mendoza Renier Editor
University of the Philippines Diliman, PHILIPPINES
Competing Interests: The authors have declared that no competing interests exist.

* E-mail: kiattisak.pr@wu.ac.th
2023
30 6 2023
18 6 e028755610 2 2023
7 6 2023
© 2023 Prathom, Jampeepan
2023
Prathom, Jampeepan
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Compartment models are implemented to understand the dynamic of a system. To analyze the models, a numerical tool is required. This manuscript presents an alternative numerical tool for the SIR and SEIR models. The same idea could be applied to other compartment models. The result starts with transforming the SIR model to an equivalent differential equation. The Dirichlet series satisfying the differential equation leads to an alternative numerical method to obtain the model’s solutions. The derived Dirichlet solution not only matches the numerical solution obtained by the fourth-order Runge-Kutta method (RK-4), but it also carries the long-run behavior of the system. The SIR solutions obtained by the RK-4 method, an approximated analytical solution, and the Dirichlet series approximants are graphically compared. The Dirichlet series approximants order 15 and the RK-4 method are almost perfectly matched with the mean square error less than 2 × 10−5. A specific Dirichlet series is considered in the case of the SEIR model. The process to obtain a numerical solution is done in the similar way. The graphical comparisons of the solutions achieved by the Dirichlet series approximants order 20 and the RK-4 method show that both methods produce almost the same solution. The mean square errors of the Dirichlet series approximants order 20 in this case are less than 1.2 × 10−4.

The author(s) received no specific funding for this work. Data AvailabilityAll relevant data are within the manuscript and its Supporting information files.
Data Availability

All relevant data are within the manuscript and its Supporting information files.
==== Body
pmcIntroduction

The SIR model is a simple system of nonlinear differential equations that has a rich dynamic. This model is an example of compartment models mainly studied in epidemiology [1–3]. The SIR model and its adaptive versions have been used broadly to investigate the dynamic of a considered system. The models are applicable in public health to intervene and control a disease’s transmission [1, 3–12]. Once a compartment model appears, a tool is required to analyze the model. Understanding the flow of the model would lead the researchers to an assumption to control the disease spreading. There are several tools available to handle the compartment models depending on the situation we are dealing with. Some researchers have sometimes used a computer model to assess certain situations [13, 14]. Others have implemented a numerical method to investigate the models. Since computer models can be utilized only in some circumstances, a numerical method would be an alternative tool. Some recent studies of the compartment models have focused on the fractional-order compartment models by considering the models under the fractional-order α, where α ∈ (0, 1). As α approaches 1, the models become the classical compartment models. The fractional-order SIR model has been demonstrated to have biological significance when α ≥ 0.8 [15].

The Dirichlet series is one of the interesting well-known series in the branch of complex analysis since it is a convergent series. The general Dirichlet series is in the form ∑i=0∞aie-λit,

where 0 = λ0 < λ1 < … and ai∈C\{0}. The started term a0 is the approached value of the series as t goes to infinity. The Dirichlet series solution satisfying some specific differential equations has been studied in [16–18]. Investigating the Dirichlet series satisfying the compartment models would lead us to a new method to obtain the models’ solution which will be mainly discussed in this paper.

This manuscript considers the classical SIR and SEIR models which are equivalent to the fractional-order α when α approaches one. The existence and uniqueness of the models’ solutions are confirmed by a well-known theorem appearing in many differential equations books. However, their exact non-parametric solutions, or what we can call the closed-form formulas have been unknown. The results in [19] have illustrated the exact analytical solution of the SIR model in the parametric form. Although the parametric analytical solution has been derived, a numerical method is still required to obtain the approximate solutions S(t), I(t), R(t) at any time t. A well-known numerical method to obtain an approximate solution of any compartment model is the fourth-order Runge-Kutta method (RK-4) which is an iterative method depending on a given step size and it does not provide the long-run behavior of the solution. An approximate closed-form solution on the other hand may give the long-run behavior with no step size needed. This rises up researchers’ interest to derive an approximate closed-form solution of the SIR and its related models. Some approximate closed-form solutions of the SIR model have been studied in [20–23]. Considering the solution in the form of a power series expression is one simple way to obtain an approximate closed-form solution of a nonlinear system of differential equations. However, the power series solution as a function of time t diverges as t approaches infinity. The power series solution of the SIR model has been studied by [24]. A multistage technique repeating the same order of power series approximation with updated initial conditions is a method to get a numerical solution of the SIR model [25]. Using power series approximation or repeatedly using a fixed order of power series approximation can not eventually reach the long-run behavior of the system due to the divergence of the power series. A power series approach together with a defined approximant is presented in [20]. Padé approximant are the other two numerical methods that can lead the approximated solution to the long-run behavior [26, 27]. Considering the Dirichlet series solution would be an alternative path to overcome the limitation mentioned above.

SIR model

In the SIR model, the rates that the proportions of susceptible S, infectious I, and recovered R change over time are described by dSdt=-βIS, (1)

dIdt=βIS-γI, (2)

dRdt=γI, (3)

where the initial conditions are assumed to be S(0) = S0 > 0, I(0) = I0 > 0, R(0) = 0, and β,γ∈R0+. Another notation used in this model is S∞:=limt→∞S(t).

Note that the solution functions S, I, R are single variable real-valued functions that S,I,R:R0+→[0,1]. In addition, the total proportion at any time t is assumed to satisfy S(t) + I(t) + R(t) = 1.

Eqs (1)-(3) can be combined to be in the form of a Bernoulli differential equation presented in [19, 21]. Solving the Bernoulli differential equation leads to an approximate analytical solution of the model. Here, we present an approach to derive the approximate analytical solution without transforming Eqs (1)-(3) into the Bernoulli differential equation. A result derived from this approach will be used further.

We first rewrite (1) as I=-1β(S′S). (4)

When we substitute this equation in (2), we get I′=-(βS-γ)(1β)(S′S). (5)

Differentiating both sides of (4) with respect to t and then equating with (5), we have SS′′-S′2-(βS-γ)SS′=0. (6)

The solution of this differential equation can be found by defining the function ϕ as 1ϕ:=1SdSdt=ddtlnS. (7)

Note that our approach of defining ϕ above is different from the approach presented in [19]. After using the notation of ϕ with (6), we have 1ϕ2dϕdS=-β+γS.

By solving this first-order differential equation and then using (7), we get ddtlnS=1ϕ=β(S-S0-I0)-γln(SS0) (8)

which is equivalent to 1β(S-S0-I0)-γln(SS0)dlnS=dt. (9)

The exact solution of (9) is called the analytical solution. By using power series approximation of ln(S/S0), we obtain an approximate solution called a semi-analytical solution. Since we assumed that S0 + I0 = 1 and the term ln(SS0) can be approximated by SS0-1, it follows that (9) becomes 1elnS(β-γS0)+γ-βdlnS=dt. (10)

We note here that the assumption I0 > 0 implies S0 < 1. It follows that the terms βS0 − γ and γ − β cannot be zero at the same time. Now, we consider two cases.

Case 1: β = γ. Integrating both sides of 10 gives S(t)=S01+β(1-S0)t. (11)

Case 2: β ≠ γ By integrating both sides of 10, we simply get S(t)=(β-γ)S0β(1-S0)e(β-γ)t+βS0-γ. (12)

Once we obtain the approximate solution S(t), we can derive the corresponding solutions of I(t) (using (4) and (8)) and R(t) = 1 − S(t) − I(t) as follows: I(t)=-S(t)+1+γβln(S(t)S0), (13)

R(t)=γβln(S0S(t)). (14)

The set of (S(t), I(t), R(t)) for t ≥ 0 presented in (12)-(14) is called a semi-analytical solution.

Dirichlet series solutions of the SIR model

Let u(t) = ln(S(t)/S0). Then, (8) is equivalent to u′+γu+β=βS0eu. (15)

Since this first-order differential equation has no exact non-parametric analytical solution, we consider its Dirichlet series solution. We first set up some notations for ease of computation. Considering 0 = λ0 < λ1 < … and ai∈C\{0}, let u(t)=∑i=0∞aie-λit=∑i=0n-1aie-λit+∑i=n∞aie-λit=:u1,n(t)+u2,n(t), (16)

and let F(u,u′)=u′+γu-βS0eu+β. (17)

Then (15) is equivalent to F(u, u′) = 0. By Taylor series expansion we have, 0=F(u,u′)=F(u1,n+u2,n,u1,n′+u2,n′)=F(u1,n,u1,n′)+An, (18)

where An=u2,n(γ-βS0eu1,n)+u2,n′-(u2,n)2βS0eu1,n2+⋯=ane-λnt(γ-βS0ea0-λn)+termswithhigherexponents. (19)

The following result provides properties of the Dirichlet series solution of (15) which is equivalent to (18).

Theorem 1 Let ln(S(t)/S0)=u(t)=∑i=0∞aie-λit , 0 = λ0 < λ1 < …, ai∈C\{0} be a Dirichlet series satisfying (15). Then,

1 − S∞ + (γ/β)ln(S∞/S0) = 0,

λ1 = γ − βS∞,

a0 = ln(S∞/S0) and a1=-a0-∑i=2∞ai,

for n ≥ 2 if 0≠ln=∑j=2n∑i1+⋯+ij=ni1,…,ij>0ai1⋯aij/j!, then λn = nλ1 and an=-βS∞ln(n-1)λ1.

Proof. It is obvious that a0 = ln(S∞/S0) and u(0)=0=∑i=0∞ai. Consider (18) with n = 1, we see that F(u1,1,u1,1′)=γa0-βS0ea0+β is a constant and the first term of A1 is non-constant, a1e-λ1t(γ-βS0ea0-λ1). This implies that 0=γa0-βS0ea0+β and λ1=γ-βS0ea0=γ-βS∞. The former is equivalent to 1-S∞+(γ/β)ln(S∞/S0)=0. (20)

Consider (18) with n = 2, observe that the first term (least exponent) of A2 is a2e-λ2t(γ-βS0ea0-λ2)≠0 since a2 ≠ 0 and λ1 < λ2. That means that this term would be canceled by a term in F(u1,2,u1,2′). Note that the exponents of the terms in F(u1,2,u1,2′) are linear combination of λ1. This implies that λ2 = nλ1 for some integer n ≥ 2. Since F(u1,2,u1,2′) has the term -βS∞a122e-2λ1t=-βS∞l2e-2λ1t≠0, it follows that λ2 = 2λ1 and then the sum of coefficients of the terms e-λ2t and e-2λ1t must be zero, i.e., a2=-βS∞l2λ1. Now consider (18) with a fixed n > 2 and assume that λk = kλ1 for 1 ≤ k ≤ n − 1. The first term of An is ane-λnt(γ-βS0ea0-λn)≠0 since an ≠ 0 and λ1 < λn. This means that the exponent −λnt must appear in F(u1,n,u1,n′). Note that F(u1,n,u1,n′)=-∑i=0n-1iλ1aie-iλ1t+γ∑i=0n-1aie-iλ1t-βS∞e∑i=1n-1aie-iλ1t+β. (21)

Since the exponent −λnt must be equal to −jλ1t for some j ≥ n, it implies that the term with the exponent −λnt can be extracted from the term -βS∞e∑i=1n-1aie-iλ1t, see (21). By the series expansion, the term -βS∞e∑i=1n-1aie-iλ1t can be written as -βS∞(1+∑i=1n-1aie-iλ1t+(∑i=1n-1aie-iλ1t)22!+(∑i=1n-1aie-iλ1t)33!+⋯)

which contains the term of exponent nλ1 in the form -βS∞(∑j=2n∑i1+⋯+ij=ni1,…,ij>0ai1⋯aij/j!)e-nλ1t=-βS∞lne-nλ1t≠0.

Hence, λn = nλ1 and also an=-βS∞ln(n-1)λ1. We have done the proof.

Note that Eqs (13), (14) and (20) imply that I∞=limt→∞I(t)=0andR∞=limt→∞R(t)=1-S∞.

This means that the solution set of the SIR model, {(S(t),I(t),R(t))|S(0)=S0,I(0)=I0,R(0)=0},

approaches (S∞, 0, 1 − S∞). Consider a Dirichlet series approximant which is an approximate form of the solution u(t), say uN(t)=∑i=0Naie-λit

for some N∈N. By the previous theorem, to obtain an we need to get ln which depends on a1, …, an−1 and a1=a0-∑i=2nai. It implies that li for each i = 3, …, N may not be unique and the values of a1, …, an can be non-unique complex numbers. This provides many approximate forms of the complex Dirichlet series u(t) which are not applicable to computer programming. We need the series u(t) to be a series of real coefficients. The results in Theorem 1 induce us to consider a specific Dirichlet series solution whose exponents satisfy λn = nλ1 for n∈N. The next results are more applicable for achieving the real Dirichlet series u(t).

Lemma 2 If f(t)=e∑i=1∞aie-λit and ∑i=0∞ai=0 , then the nth derivative of f with respect to t at t = 0 satisfies f(n)(0)=(-1)ne-a0∑k=1nk−1n−1cn-kbk,

where ck=k−10ck−1b1+k−11ck−2b2+⋯+k−1k−1c0bk for k ≥ 1, c0 = 1 and bk=∑i=1∞aiλik.

The proof simply follows from the general Leibniz rule and the mathematical induction.

Theorem 3 Let ln(S(t)/S0)=u(t)=∑i=0∞aie-λit , 0 = λ0 < λ1 < …, ai∈R\{0} be a Dirichlet series satisfying (15). Let bn=∑i=1∞aiλin for n ≥ 0. Then, b0=-a0,b1=βI0,bn=γbn-1-βS0∑i=1n-1i−1n−2cn-1-ibi

for n = 2, …, where ck=k−10ck−1b1+⋯+k−1k−1c0bk for k ≥ 1 and c0 = 1.

Proof. By the result in Theorem 1, we have b0 = −a0. By plugging the Dirichlet series u(t) into (15), we get β+γa0+∑i=1∞ai(γ-λi)e-λit=βS∞e∑i=1∞aie-λit. (22)

Eq (22) at t = 0 implies that b1 = βI0. By taking the (n − 1)th derivative both sides of (22) with respect to t and then plugging t = 0 together with the result of Lemma 2, we obtain that ∑i=1∞(-1)n-1aiλin-1(γ-λi)=βS∞(-1)n-1e-a0∑i=1n-1n-2i-1cn-1-ibi.

This implies the result of this theorem.

Dirichlet series approximants of the basic SIR model

Followed by Theorem 1, it is intuitive to consider the case of λi = iλ1. This implies by setting in Theorem 3 that ∑i=1∞ai=-a0,∑i=1∞iai=b1/λ1,⋮∑i=1∞inai=bn/λ1n,⋮

which can be written in the matrix form as follows: [11⋯1⋯12⋯n⋯1222⋯n2⋯⋮⋮⋮⋮⋯1n2n⋯nn⋯⋮⋮⋮⋮⋮][a1a2a3⋮an⋮]=[-a0b1/λ1b2/λ12⋮bn/λ1n⋮],

where a0 = ln(S∞/S0), S∞ can be derived by (20), and b1, b2, … can be obtained by the result of Theorem 3.

Now consider an approximate form of the Dirichlet series uN(t)=∑i=0Naie-iλ1t satisfying (15). Then by the previous results, we get a0=ln(S∞/S0),1=S∞-γβln(S∞S0),λ1=γ-βS∞,

and a1, …, aN can be achieved by solving [11⋯112⋯N1222⋯N2⋮⋮⋮⋮1N-12N-1⋯NN-1][a1a2a3⋮aN]=[-a0b1/λ1b2/λ12⋮bN/λ1N-1]. (23)

The existence of the solution a1, a2, …, aN is due to the existence of the inverse of the Vandermonde matrix [28]. Note here that this matrix equation is solvable when N, a0, λ1, b1, …, bN are known. After uN(t) is obtained by solving the matrix equation above, an approximation of the solution S(t), say SN(t), is derived by SN(t)=S0euN(t). (24)

The other approximations of the solutions I(t) and R(t), say IN(t) and RN(t) respectively, are as follows: IN(t)=-SN(t)+1+γβln(SN(t)S0) (25)

RN(t)=γβln(S0SN(t)). (26)

It is implied by Theorem 3 that ck can be written in terms of b1, …, bk. The process to get ck presented in Theorem 3 is so applicable that any mathematical programming can perform the task. The terms bk can be derived after the terms c0, …, ck−2 were known.

Numerical simulations for the SIR model

In this subsection, we present an algorithm to obtain the Dirichlet series approximants SN(t), IN(t), RN(t). According to the previous results, the process of achieving the approximant of the SIR model is as follows:

Step 1: Setting c0 = 1, then ck, k∈N, can be written in terms of b1, b2, …, bk by using ck=k−10ck−1b1+k−11ck−2b2+⋯+k−1k−1c0bk.

Step 2: Solving 1 − S∞ + (γ/β)ln(S∞/S0) = 0 for S∞, setting a0 = ln(S∞/S0), λ1 = r − βS∞ and b0 = −a0.

Step 3: Setting b1 = βI0, then bk can be written in terms of b1, …, bk−1 by using Steps 1 and bk=γbk-1-βS0∑i=1k-1i−1k−2ck-1-ibi. This step gives the values of b1, b2, …, bN.

Step 4: Solving the matrix Eq (23) for a1, a2, …, aN. In this step we get uN(t)=∑i=0Naie-iλ1t.

Step 5: Plotting SN(t), IN(t), RN(t) by using Eqs (24–26)

We first show manually how to get some cn by performing Step 1 c0=1,c1=b1,c2=c1b1+b2=b12+b2,c3=b13+3b1b2+b3,c4=(b13+3b1b2+b3)b3+3(b12+b2)b2+3b1b3+b4,

and then we manually derive bn by implementing Step 3 and the result in Step 1 b0=-a0=-ln(S∞/S0),b1=βI0,b2=(γ-βS0)b1=(γ-βS0)βI0,b3=γβI0(γ-βS0)-β2S0I0(γ+βI0-βS0),b4=βI0(γ-βS0)(γ2-2γβS0-3β2I0S0+β2S02)-γβ3I02S0-β4I03S0+β4S02I02.

Fig 1 shows the comparison of the SIR solutions obtained by RK-4 method, Dirichlet series approximants (S6(t), S10(t), S15(t)), and the semi-analytical solution. Even though the graph of the semi-analytical solution using (12) is close to the numerical solution at the beginning time interval, it has no biological meaning after 21 days passed. The infectious proportion is negative and the recovered fraction is higher than 1 when 21 days have passed, see Fig 1. The approximants SN(t), IN(t), RN(t) are close to the numerical solution on the beginning and ending time intervals. As shown in Fig 1, the better approximation comes with the larger N. Figs 1 and 2 also illustrate that a small positive integer N is enough to achieve a suitable numerical approximation of the SIR solution.

10.1371/journal.pone.0287556.g001 Fig 1 The solution graph comparison.

Under the assumed parameters that β = 0.5, γ = 0.2 and the initial conditions S0 = 0.99, I0 = 0.01, R0 = 0, the three figures show the graphs of the proportions of susceptible, infectious, and recovered cases. The figures show the comparison of the approximate solution obtained by the numerical method (RK-4), the Dirichlet series approximants SN(t), IN(t), RN(t), for N = 6, 10, 15, and the approximate analytical using (12).

10.1371/journal.pone.0287556.g002 Fig 2 The Dirichlet approximations order 15 Vs RK-4 method.

Under the assumed parameters related to the COVID-19 situation [2] β = 0.35, γ = 0.14, and the initial conditions S0 = 0.99, I0 = 0.01, R0 = 0, the figure compares the graphs of the SIR model solutions obtained by the numerical method (RK-4) in the solid lines and the Dirichlet series approximants S15(t), I15(t), R15(t) in the dashed lines.

Here we show u15(t) corresponding to the given parameters and initial conditions in Fig 2. u15(t)=-2.23526+0.667927e-1.54406t-10.1641e-1.44112t+72.3135e-1.33818t-319.195e-1.23524t+978.052e-1.13231t-2205.03e-1.02937t+3781.85e-0.926433t-5030.29e-0.823496t+5240.03e-0.720559t-4284.69e-0.617622t+2736.93e-0.514685t-1348.28e-0.411748t+500.112e-0.308811t-133.847e-0.205874t+23.7809e-0.102937t

By using this real Dirichlet series u15(t) with (24)-(26), we get the approximated solution of the SIR model. The results above also give the long-run behavior of the considered SIR model in Fig 2 as follows: (S∞,I∞,R∞)=(S∞,0,1-S∞)=(0.105894,0,0.894106).

As seen in Figs 1 and 2, the approximations S15(t), I15(t), R15(t) are perfectly matched to the numerical solution obtained by the RK-4 method.

We note here that the approximate closed-form derived by our method presents an approximate long-run behavior of the system while the iterative RK-4 method and the exact parametric solution in [19] do not present it.

Error comparison of the SIR model

Since the exact values of the solutions S(t), I(t), and R(t) at any time t are unknown because the model has no closed-form solution, we illustrate the error of our method comparing with the RK-4 method. Fig 3 shows the errors of S15(t), I15(t), R15(t) compared with the RK-4 method. The graphs illustrate the gaps between the solutions from the RK-4 method and our method; i.e., the graphs show the differences S15(t) − S(t), I15(t) − I(t), and R15(t) − R(t) at any time t. The mean squared errors of S15, I15(t), R15(t) compared with the numerical solutions S(t), I(t), R(t) are 5.88606 × 10−6, 5.81967 × 10−6 and 0.0000187405, respectively.

10.1371/journal.pone.0287556.g003 Fig 3 Error of the solutions S15(t), I15(t), R15(t) compared with the solutions obtained from the RK-4 method.

A special Dirichlet series solutions of the SEIR model

In this section, we consider another compartment model consisting of the susceptible population S, the exposed population E, the infectious population I, and the recovered population R called the SEIR model whose system of differential equations is as follows: dSdt=-βIS, (27)

dEdt=βIS-σE, (28)

dIdt=σE-γI, (29)

dRdt=γI, (30)

where the initial conditions are assumed to be S(0) = S0 > 0, E(0) = E0 > 0, I(0) = I0 > 0, R(0) = 0, and β,σ,γ∈R0+. Additionally, the total proportion at any time t is assumed to satisfy S(t) + E(t) + I(t) + R(t) = 1.

Let us say that we want an approximated solution R(t) of this model to be a specific Dirichlet series in the form: R(t)=∑i=0∞hie-iλ1t,hi,λ1∈R\{0}. (31)

By plugging (31) in (30) and (29), it is easy to see that the solutions I(t) and E(t) can be written in a similar Dirichlet form of R(t) and the constant terms of the Dirichlet expansion of I(t) and E(t) are zero. This implies that I∞:=limt→∞I(t)=0andE∞:=limt→∞E(t)=0.

We also simply obtain S∞+R∞:=limt→∞S(t)+limt→∞R(t)=1. (32)

By substituting (30) in (27), we get dS=-βSdR.

Solving this separable differential equation leads to ln(S(t)S0)=-βγR(t). (33)

By (32) and (33), we have h0=R∞=-γβln(1-R∞S0). (34)

Eq (31) with t = 0 implies that ∑i=1∞hi=-h0.

Plugging (31) in (30) and setting t = 0 lead to ∑i=1∞ihi=-γI0λ1.

Now setting Dn ≔ γ(σE(n)(0) − γI(n)(0)), where E(n)(0) and I(n)(0) are the nth derivative of E(t) and I(t) at t = 0, respectively, for n∈N∪{0}. Taking the first derivative of (30) with respect to t and then plugging t = 0 lead us to ∑i=1∞i2hi=γ(σE0-γI0)λ12=D0λ12.

If we keep taking derivatives till the nth derivative of R(t) with respect to t and plugging t = 0, we should get ∑i=1∞inhi=(-1)nγ(σE(n-2)(0)-γI(n-2)(0))λ1n=(-1)nDn-2λ1n.

An approximated solution R(t) could be considered in the form RN(t)=∑i=0Nhie-iλ1t,

where h0 is obtained by (34) and other hi are achieved by solving the matrix equation: [11⋯112⋯N1222⋯N2⋮⋮⋮⋮1N-12N-1⋯NN-1][h1h2h3⋮hN]=[-h0-γI0/λ1D0/λ12⋮(-1)N-1DN-3/λ1N-1]. (35)

Note that Dn can be recursively solved by Dn=γ(σE(n)(0)-γI(n)(0)), (36)

I(n)(0)=σE(n-1)(0)-γI(n-1)(0), (37)

E(n)(0)=β(IS)(n-1)(0)-σE(n-1)(0), (38)

S(n)(0)=-β(IS)(n-1)(0), (39)

(IS)(n-1)(0)=∑k=0n-1(n−1k)I(n-1-k)(0)S(k)(0). (40)

Eqs (35)-(40) give us an approximated solution RN(t, λ1). After the solution RN(t,λ1)=∑i=0Nhie-iλ1t

is derived, other solutions can be achieved as follows: IN(t,λ1)=1γddtRN(t,λ1), (41)

EN(t,λ1)=1σ(γIN(t,λ1)+ddtIN(t,λ1)), (42)

SN(t,λ1)=1-EN(t,λ1)-IN(t,λ1)-RN(t,λ1). (43)

The approximated solution (SN(t), EN(t), IN(t), RN(t)) of the SEIR is obtained by choosing a value of λ1.

Numerical simulations for the SEIR model

In this subsection, we present a method to achieve an approximated solution of the SEIR model by using the results of the previous subsection. Given the system of differential Eqs (27)-(30) with known parameters β, σ, γ and initial conditions S0, E0, I0, R0, the process to derive the approximated solution is the following:

Step 1: Derive Dn by using (36)-(40).

Step 2: Solving Eq (34) for h0,

Step 3: Solving the matrix Eq (35) for h1, h2, …, hN. In this step we get RN(t,λ1)=∑i=0Nhie-iλ1t.

Step 4: Fixed λ1 and Plotting SN(t), EN(t), IN(t) and RN(t) by using Eqs (41–43).

Note that each Dn is the function of know parameters β, σ, γ, S0, E0, I0. Here are some examples of Dn expressions: D0=γ(σE0-γI0),D1=γσ(βI0S0-σE0)-γ2(σE0-γI0)D2=(γ3+γσβS0)(σE0-γI0)-(γ2σ+γσ2)(βI0S0-σE0)-γσS0(βI0)2.

To illustrate the results generated by this process, we consider the SEIR model (27)-(30) with the initial conditions S0 = 0.98, E0 = 0.01, I0 = 0.01, and R0 = 0 under the assumed parameters obtained by [29], β = 0.75 per day, σ = 1/3 per day, γ = 1/8 per day. The comparisons are demonstrated by Figs 4 and 5. Note that in those figures, the value of λ1 is chosen to be 0.1. In Fig 4, the graphs of the Dirichlet series approximants SN(t), EN(t), IN(t), and RN(t) for N = 4, 12, 20 created by our algorithm are compared with the solution from the RK-4 method. As we can see in the figure, the Dirichlet series approximants get closer to the numerical solution obtained by the RK-4 method when N gets larger. Fig 5 is shown to conclude that the Dirichlet series approximants order 20, S20(t), E20(t), I20(t) and R20(t) are almost perfectly match the numerical method RK-4. The solution R20(t) obtained by the algorithm is as follows: R20(t)=0.997535-0.906981e-2t+18.0889e-1.9t-171.28e-1.8t+1023.69e-1.7t-4330.52e-1.6t+13780.6e-1.5t-34219e-1.4t+67872.2e-1.3t-109161e-1.2t+143670e-1.1t-155430e-t+138262e-0.9t-100721e-0.8t+59540.1e-0.7t-28099.4e-0.6t+10297e-0.5t-2786.09e-0.4t+499.369e-0.3t-40.602e-0.2t-3.74085e-0.1t.

10.1371/journal.pone.0287556.g004 Fig 4 The comparison of solutions SN(t), EN(t), IN(t), RN(t) for N = 4, 12, 20 versus the solution obtained by the numerical method RK-4.

10.1371/journal.pone.0287556.g005 Fig 5 The comparison of the solution (S20, E20, I20, R20) by our method versus the numerical method RK-4.

This real Dirichlet series solution R20(t) leads to other solutions I20, E20, and S20 by using (41)-(43). Note that for each N the solution RN(t) gives the long-run behavior R∞. In the considered situation, we have R∞ = 0.997535 which is obtained by solving Eq (34). We also can note here that the constant 0.997535 is the constant term of RN(t) for any integer N. It can be concluded that under the given initial conditions and parameters in this subsection, the long-run behavior of the SEIR solution is (S∞,E∞,I∞,R∞)=(0.002465,0,0,0.997535).

Error comparison of the SEIR model

The accuracy of our method can be observed by investigating the mean square error. Fig 6 shows the errors of S20(t), E20(t), I20(t), R20(t) compared with the numerical solution obtained by the RK-4 method. The figure shows the graphs of the differences S15(t) − S(t), E15(t) − E(t), I15(t) − I(t), and R15(t) − R(t) at any time t. The mean square errors of S20(t), E20(t), I20(t), R20(t) compared with the numerical solutions S(t), E(t), I(t), R(t) are 0.000111401, 0.0000150904, 0.0000287791 and 0.0000675307, respectively.

10.1371/journal.pone.0287556.g006 Fig 6 Error of the solutions S20(t), E20(t), I20(t), R20(t) compared with the solutions S(t), E(t), I(t), R(t) obtained from the RK-4 method.

Conclusions

The presented method needs no step-size interpolation and also shows the long-run behavior of the approximate solutions; however, the exact parametric solution (SIR case) and the solution from the RK-4 method (both SIR and SEIR cases) cannot give the long-run behavior of the system since their processes need step-size interpolation. The exact parametric analytical solution (SIR case) presented in [19] does not give the closed-form solution as a function of time and it still needs a numerical method to get the SIR solution at any time t. The approximate closed-form Dirichlet series solution can be an alternative direct method to get the approximate solution of the SIR and SEIR models.

The theoretical results in this paper provide a numerical algorithm in order to obtain a numerical solution of the SIR and SEIR models. The derived Dirichlet series approximants obviously preserve the positive property of the SIR and SEIR solutions. The comparisons of the solution graphs shown in Figs 1, 2, 4, and 5 assure that the presented method in this manuscript can be used as an alternative method. The SIR solution in the Dirichlet series form obviously gives the long-run behavior (S∞, 0, 1 − S∞) which depends on the given initial conditions and the model’s parameters. As a consequence of the approach to the SIR model, a similar process is used with the SEIR model. The assumed characteristic of the solution R(t) as a real Dirichlet series leads us to the approximated solutions of the SEIR model. The solutions SN(t), EN(t), IN(t), RN(t)) obtained by our method is almost perfectly matched the solution from the RK-4 method. The solutions also obviously preserve the positive property of the SEIR model and carry the long-run behavior (S∞, E∞, I∞, R∞) which depends on the initial conditions and parameters. The method presented in this manuscript is not only useful for obtaining the SIR and SEIR models but it also can be used in many kinds of nonlinear dynamical systems.

Supporting information

S1 File Wolfram mathematica codes.

There are no raw data used in this paper. All figures in this manuscript are generated by Wolfram Mathematica posted at https://github.com/Kiattisak-Prathom/Mathematica-Codes-for-PLOS-ONE.

(ZIP)

Click here for additional data file.

The authors appreciate the journal’s anonymous reviewers and editors for their insightful criticism and recommendations.
==== Refs
References

1 Weiss HH . The SIR model and the foundations of public health. Materials matematics. 2013:0001–17.
2 Cooper I , Mondal A , Antonopoulos CG . A SIR model assumption for the spread of COVID-19 in different communities. Chaos, Solitons & Fractals. 2020 Oct 1;139 :110057. doi: 10.1016/j.chaos.2020.110057 32834610
3 Acemoglu D , Chernozhukov V , Werning I , Whinston MD . Optimal targeted lockdowns in a multigroup SIR model. American Economic Review: Insights. 2021 Dec;3 (4 ):487–502.
4 Wintachai P , Prathom K . Stability analysis of SEIR model related to efficiency of vaccines for COVID-19 situation. Heliyon. 2021 Apr 1;7 (4 ):e06812. doi: 10.1016/j.heliyon.2021.e06812 33880423
5 Biswas MH , Paiva LT , de Pinho MD . A SEIR model for control of infectious diseases with constraints. Mathematical Biosciences and Engineering. 2014 Aug 1;11 (4):761–84. doi: 10.3934/mbe.2014.11.761
6 Mwalili S , Kimathi M , Ojiambo V , Gathungu D , Mbogo R . SEIR model for COVID-19 dynamics incorporating the environment and social distancing. BMC Research Notes. 2020 Dec;13 (1 ):1–5. doi: 10.1186/s13104-020-05192-1 31898526
7 López L , Rodo X . A modified SEIR model to predict the COVID-19 outbreak in Spain and Italy: simulating control scenarios and multi-scale epidemics. Results in Physics. 2021 Feb 1;21 :103746. doi: 10.1016/j.rinp.2020.103746 33391984
8 Ragusa MA, Razani A. Weak solutions for a system of quasilinear elliptic equations. arXiv preprint arXiv:2006.05262. 2020 Jun 6.
9 Goodarzi Z , Mokhtarzadeh MR , Pournaki MR , Razani A . A Note On Periodic Solutions Of Matrix Riccati Differential Equations. Applied Mathematics E-Notes. 2021;21 :179–86.
10 Razani A . An existence theorem for ordinary differential equation in Menger probabilistic metric space. Miskolc Mathematical Notes. 2014;15 (2 ):711–6. doi: 10.18514/MMN.2014.640
11 Bjørnstad ON , Shea K , Krzywinski M , Altman N . The SEIRS model for infectious disease dynamics. Nature methods. 2020 Jun 1;17 (6 ):557–9. doi: 10.1038/s41592-020-0856-2 32499633
12 Trejo I , Hengartner NW . A modified Susceptible-Infected-Recovered model for observed under-reported incidence data. PloS one. 2022 Feb 9;17 (2 ):e0263047. doi: 10.1371/journal.pone.0263047 35139110
13 Cuspilici A , Monforte P , Ragusa MA . Study of Saharan dust influence on PM10 measures in Sicily from 2013 to 2015. Ecological Indicators. 2017 May 1;76 :297–303. doi: 10.1016/j.ecolind.2017.01.016
14 Duro A, Piccione V, Ragusa MA, Veneziano V. New environmentally sensitive patch index-ESPI-for MEDALUS protocol. InAIP Conference Proceedings 2014 Dec 10 (Vol. 1637, No. 1, pp. 305-312). American Institute of Physics.
15 Salman SM . On a discretized fractional-order SIR model for Influenza A viruses. Prog. Fract. Differ. Appl. 2017;3 (2 ):163–73. doi: 10.18576/pfda/030207
16 Das AG , Lahiri BK . Dirichlet series solutions of differential equations. Rendiconti del Circolo Matematico di Palermo. 1984 Sep;33 :425–35. doi: 10.1007/BF02844504
17 Sachdev PL , Bujurke NM , Pai NP . Dirichlet series solution of equations arising in boundary layer theory. Mathematical and computer modelling. 2000 Nov 1;32 (9 ):971–80. doi: 10.1016/S0895-7177(00)00183-7
18 Laohakosol V , Phuksuwan O , Prathom K . Dirichlet series solutions of generalized Riccati equations. European Journal of Mathematics. 2015 Mar;1:170–85. doi: 10.1007/s40879-014-0021-5
19 Harko T , Lobo FS , Mak M . Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates. Applied Mathematics and Computation. 2014 Jun 1;236 :184–94. doi: 10.1016/j.amc.2014.03.030
20 Barlow NS , Weinstein SJ . Accurate closed-form solution of the SIR epidemic model. Physica D: Nonlinear Phenomena. 2020 Jul 1;408 :132540. doi: 10.1016/j.physd.2020.132540 32362697
21 Heng K , Althaus CL . The approximately universal shapes of epidemic curves in the Susceptible-Exposed-Infectious-Recovered (SEIR) model. Scientific Reports. 2020 Nov 9;10 (1 ):19365. doi: 10.1038/s41598-020-76563-8 33168932
22 Yıldırım A , Cherruault Y . Analytical approximate solution of a SIR epidemic model with constant vaccination strategy by homotopy perturbation method. Kybernetes. 2009 Oct 16;38 (9 ):1566–75. doi: 10.1108/03684920910991540
23 Izadi M , Seifaddini M , Afshar M . Approximate solutions of a SIR epidemiological model of computer viruses. Advanced Studies: Euro-Tbilisi Mathematical Journal. 2021 Dec;14 (4 ):203–19.
24 Srivastava HM , Area Carracedo IC , Nieto JJ . Power-series solution of compartmental epidemiological models. Mathematical Biosciences and Engineering. 2021 Apr 12. doi: 10.3934/mbe.2021163 34198385
25 Akindeinde SO . A new multistage technique for approximate analytical solution of nonlinear differential equations. Heliyon. 2020 Oct 1;6 (10 ):e05188. doi: 10.1016/j.heliyon.2020.e05188 33088955
26 Vazquez-Leal H , Benhammouda B , Filobello-Nino U , Sarmiento-Reyes A , Jimenez-Fernandez VM , Garcia-Gervacio JL , et al . Direct application of Padé approximant for solving nonlinear differential equations. SpringerPlus. 2014 Dec;3 (1 ):1–1. doi: 10.1186/2193-1801-3-563 24422185
27 Baker GA , Baker GA Jr , Graves-Morris P , Baker SS . Pade Approximants: Encyclopedia of Mathematics and It’s Applications, Vol. 59 Baker George A. Jr. , Peter Graves-Morris . Cambridge University Press; 1996 Jan 26.
28 Turner LR . Inverse of the Vandermonde matrix with applications. 1966 Aug 1.
29 Carcione JM , Santos JE , Bagaini C , Ba J . A simulation of a COVID-19 epidemic based on a deterministic SEIR model. Frontiers in public health. 2020 May 28;8 :230. doi: 10.3389/fpubh.2020.00230 32574303