Numerical Solution of an Interval-Based Uncertain SIR (Susceptible–Infected–Recovered) Epidemic Model by Homotopy Analysis Method

: This work proposes an interval-based uncertain Susceptible–Infected–Recovered (SIR) epidemic model. The interval model has been numerically solved by the homotopy analysis method (HAM). The SIR epidemic model is proposed and solved under different uncertain intervals by the HAM to obtain the numerical solution of the model. Furthermore, the SIR ODE model was transformed into a stochastic differential equation (SDE) model and the results of the stochastic and deterministic models were compared using numerical simulations. The results obtained were compared with the numerical solution and found to be in good agreement. Finally, various simulations were done to discuss the solution.


Introduction
Interval analysis is a method developed by mathematicians in the 1950s as a way of handling bounds or rounding errors and measurement errors in mathematical computation. It is useful in formulating numerical methods that yield desirable results. In short, it defines each value as a range of possibilities. This work aims to formulate interval arithmetic that solves upper and lower endpoints for the range of values of a particular function in one or more variables. These limitations are not necessarily the supremum or infimum since the exact solution of those values can be very intractable or even impossible. The treatments of interval arithmetic for real intervals of quantities with the form [u, v] = {x ∈ R : u ≤ x ≤ v}, where u = −∞ and v = ∞, are permitted. The permission is based on the fact that if one of the real intervals is infinite, we would have an unbounded interval, and if both are infinite, we would have the extended real number system. Considering the classical calculation with real numbers, simple arithmetic operations and functions on elementary intervals must initially be defined. It is after this that complicated functions can be evaluated from the basic elements. In interval arithmetic, we state the range of possible outcomes explicitly.
Thus, the results are no longer stated as numbers but as intervals, which denotes imprecise values. With the size of the intervals, we express the extent of uncertainties, which are similar to error bars to a metric. The evaluations of the outer bounds of intervals are enabled by simple arithmetic operations, for example, basic arithmetic and trigonometric. Interval arithmetic was introduced by [1] as an approach to bound rounding errors in mathematical computation. The theory of interval analysis emerged considering the computation of both the exact solution and the error term as a single entity, that is, the interval. Though a simple idea, it is a very powerful technique with numerous applications in mathematics, computer science, and engineering.
In their survey, they discussed the basic concepts of interval arithmetic and some of its extensions. They also reviewed the successful application of this theory in computer science, in particular. The authors of [2] investigated the solution of linear and nonlinear ordinary differential equations with the fuzzy initial condition. They proposed two Euler-type methods to obtain a numerical solution to the problem. They also compared their solution with existing results. They observed that the results obtained were tighter than the results from the existing method. The authors of [3] also investigated the numerical solution of n-th order fuzzy differential equations in the fuzzy environment using a homotopy perturbation method (HPM). They used triangular fuzzy convex normalized sets for the fuzzy parameter and variables.
They also compared their results obtained with the existing solution in terms of plots to show the efficiency of their method. The authors of [4] gave an overview of applications of interval arithmetic and discussed the verification methods for linear and nonlinear systems of equations. They also then discussed item software in the field and gave some historical remarks. The authors of [5] provided algorithms for computing the operations of interval arithmetic. They generated data that are sufficiently detailed to convert directly to a program to efficiently implement the interval operations. Finally, they extended these results to the case of general intervals, which are defined as connected sets of rules that are not necessarily closed. For this present work, we considered an interval-based uncertain epidemic model. A related mathematical model was proposed first for SIR transmission dynamics and then the HAM was applied to find the solution. This method employs the concept of the homotopy from topology to generate a convergent series solution of nonlinear systems. The convergent series solution of nonlinear systems was enabled by utilizing a homotopy-MacLaurin series to deal with the nonlinearity in the system. The HAM is much better than most of the existing analytic approximation method because most of the existing methods are valid only for weakly nonlinear problems [6]. It overcomes the restrictions of all other analytic approximation methods and is valid for highly nonlinear problems [6]. The HAM is always valid even if small physical parameters exist or not, it provides an easy way to guarantee the convergence of approximation series, and lastly, it provides sufficient freedom to choose the equation type of sub-problems and the base function of solutions [6]. The strength of the HAM to naturally exhibit convergence of the series solution is strange in most analytic and semi-analytic approaches to nonlinear PDEs [7]. Recently, [8] used the HAM approach to solve the SIS and SIR models of [9]. The authors of [10], extended the work of [8] to solve the SIR epidemic model in the presence of a constant vaccination strategy. The authors of [7] also applied the HAM to solve the SIR epidemic model. They obtained an explicit analytic solution of the coupled nonlinear differential equations describing the epidemic model proposed. They also compared the numerical results, which showed that the two results are in good agreement. The authors of [11] studied a new approach for solving the SIR epidemic model using the HAM that was based on dividing the entire domain into subintervals.
Other works on the homotopy analysis method with the SIR model can be found in [12][13][14][15][16][17][18][19][20]. The aim of this work was to obtain the numerical solution of an interval-based uncertain SIR epidemic model using the HAM and comparing their stochastic version. The homotopy analysis method (HAM) has been applied here to study the solution of the epidemic model under uncertain intervals. The results obtained by the HAM were compared with the approximate solution and were found to be in strong agreement. We have also developed the stochastic version of the SIR epidemic model presented in this paper in order to measure the effect of randomness of the variables in the model. To the best of our knowledge, no work has been done in the area of an interval-based uncertain SIR epidemic model and very few works have been done on the stochastic model of SIR epidemic models so far. The paper is organized as follows: in Section 2, preliminaries and basic definitions are presented. In Section 3, the presentation of the proposed model is made. In Section 4, we describe the interval-based uncertain model. In Section 5, we present the homotopy analysis approach to a non-linear system, while in Section 6, we present the solution of the SIR epidemic model by the HAM. In Section 7, we present the solution, numerical results, and a discussion on the interval-based uncertain SIR epidemic model. Section 8 showcases the stochastic version of the model. In Section 9, the graphical illustrations of our results are discussed. In Section 10, the numerical solution of the SDE model are discussed. In Section 11, we present the discussion, conclusion, and possible extensions, and finally, the references are presented.

Preliminaries
In this section, we present some notations, definitions, and preliminaries that are used further in this paper. A.
Interval Arithmetic [1] Interval arithmetic is defined on the sets of intervals, instead of sets of real numbers. Interval arithmetic defines a set of operations on intervals, as follows: where u and v are intervals.
B. Closed Interval [1] A closed interval, denoted by [m, n], is the set of real numbers given by C. Endpoint notation, interval equality [1] Two intervals, A and B, are said to be equal if they are the same sets. Hence, operationally this occurs if their corresponding endpoints are equal; A = B if A = B and A = B.
Here A, represents the left endpoint of an interval A while A represents the right endpoint of an interval A, such that A = A, A . D. Midpoint of A [1] The midpoint of A is given by E. Interval Arithmetic and Operations [1] The key point in the definition of arithmetic operations is that computing intervals are computing with sets. Let Then, the following properties hold: (i) The sum of two intervals, A and B, is the set (ii) The difference of two intervals, A and B, is the set (iii) The product of A and B is given by (iv) The quotient of A/B is defined as

Model Formulation
A population comprising three kinds of individuals, denoted by S (susceptible human), I (infected human), and R (recovered human), are considered. The susceptible human ((S(t)) is the number of susceptible humans at time t, that is, humans who are vulnerable and are yet to contract the disease but have a probability of contracting it. The infected human ((I(t)) is the population of the infected and infectious persons who have the disease and can transmit it to others, while the recovered human ((R(t)) is the population of recovered humans who cannot get the disease or transmit it, because they have natural immunity, they have recovered from the disease and are immune to re-infection, they have been placed in isolation, or they have died.
The population of susceptible humans is generated through the reduction of the rate of transmission β with the infected, such that the rate of change of the population of susceptible humans is given by the following: The rate of change of the population of infected humans is increased by the rate of transmission β with the susceptible, and reduced by the rate at which the infected population becomes isolate or recovered γ. Hence it is given by The population of recovered humans is generated by the rate at which the infected population becomes isolated or recovered. Hence it is given by Hence, the governing equation by [9] related to the present model is given by Subject to the initial conditions,

Interval-Based Uncertain Model
As mentioned in the introduction, if we assumed that the parameters involved in a model are given in terms of an interval then it will become an interval-based model and the solution has to be handled carefully. As such, let us suppose that we have the rate of transmission β and the rate at which the infected population become isolated or recovered γ in terms of intervals β = β, β and γ = γ, γ then the corresponding interval model may be written as dS dt = − βSI, with the initial conditions where S, I, and R are all now in interval form. It may be noted from the open literature that the involved parameters, such as β and γ, are usually given in term of some ranges, so we have investigated the problem considering those ranges in terms of intervals. Hence, the intervals of β and γ are taken as the following: Next, the above interval model has been solved by the homotopy analysis method (HAM). We provide, in the next section, some mathematical results.

Mathematical Results
We assume here that all parameters in Equation (5) are positive intervals. For the SIR model (5) to be meaningful biologically, we need to prove that all its stated variables are non-negative (except S) for all time, that is, the solutions of the Equation (5) with non-negative initial data will remain non-negative for all time t > 0.
Proof. We established that Equation (5) can be rewritten in the following form: The functions f 1 and f 2 are C ∞ . Thus, according to the Cauchy-Lipschitz theorem [21], Equation (5) has a unique positive solution (S, I) on the infinite time horizon, whenever S 0 > 0, I 0 > 0.  Proof. Equation (5) can be rewritten as follows: Applying the assumption that Thus, the field remains on the domain Ω.
In contrast, we show that by setting S(t) and I(t) as continuous intervals, such that S(0) = S 0 > 0 and I(0) = I 0 > 0, if I t < 0, then by the intermediate value theorem, there exists τ 1 ∈ 0, t , such that I(τ 1 ) = 0. By applying the second equation of Equation (5), we obtain I(t) = I(τ 1 )e g = 0 for t ≥ t 0 , where g is the base of − βI(t) − γI. Therefore, I(t) = 0 for t ≥ τ 1 which is a contradiction. We apply the same arguments to S(t). We show this according to Proposition 4. (5), then S(t) ≥ 0,

Proposition 4 ([22]). Suppose S(t), I(t), R(t) is a solution of Equation
If we add the first two equations of Equation (5) Therefore, Hence, there exists a point S ∞ , uniquely, 0 < S ∞ < S 0 such that I(S ∞ ) = 0 and I(S) > 0 for S ∞ < S ≤ S 0 . The point (S ∞ , 0) is called the equilibrium point of the first two Equations of (5) since both dS/dt and dI/dt vanish at t = 0. We show this according to Proposition 5.
By dividing dS dt by dR dt , which yields So that We show this according to Lemma 1.

Lemma 1 ([22]). Suppose (S(t), I(t), R(t)) be a solution of Equation (5) in the domain
Recall that from the first equation of Equation (5) and Proposition 5, we have dS dt = − βSI ≤ 0 and we say S(t) is a decreasing function, then lim t→∞ S(t) = S ∞ , such that S ∞ is a finite number. Recall also from Equation (5), the third equation dR dt = γI ≥ 0 and we say R(t) is an increasing function. Hence, by Proposition 5, lim t→∞ R(t) = R ∞ , then R ∞ is a finite number. We show this according to Lemma 2.

Lemma 2 ([22]
). If (S(t), I(t), R(t)) is a solution of Equation (5), then S(t) → S ∞ as R(t) → R ∞ as t → ∞ , such that S ∞ and R ∞ are finite numbers. Alternately, we integrate the first equation of Equation (5): We hereby present below the procedure for the HAM for the benefit of finding the numerical solution of our interval-based uncertain model. Consider a nonlinear equation where A is a linear operator, t denotes the time, and v(t) is an unknown function. Let v 0 (t) denote an initial approximation of v(t) and Z denote an auxiliary linear operator [21]. We construct the zero-order deformation equation where q ∈ [0, 1] is the embedding parameter and h = 0 is a non-zero auxiliary function. When q = 0 and q = 1, the zero-order deformation equation becomes, respectively, and ϕ(t; 1) = ϑ 0 (t).
Thus, as q increases from 0 to 1, the solution ϕ(t; q) varies continuously from the initial approximation ϑ 0 (t) of the exact solution ϑ(t). Such a kind of continuous variation is called deformation in topology. Expanding ϕ(t; p) by the Taylor series in the power series of q, we have where is the deformation derivative. If the auxiliary linear operator A, the initial approximation v 0 (t), the auxiliary parameter h I and the auxiliary function H(t) are properly chosen so that (i) the solution ϕ(t; q) of the zero-order deformation Equation (6) exists for all q ∈ [0, 1], (ii) the deformation derivative (11) exists for all m = 1, 2, . . ., (iii) the series (10) converge at q = 1, then, we have the series solution: Define the vector as → ϑ m (t) = {ϑ 0 (t), ϑ 1 (t), . . . , ϑ m (t)}. (13) According to the definition (10), the governing equation can be derived from the zeroorder deformation Equation (6). Differentiating (6) m-times with respect to the embedding parameter q, then by setting q = 0, and finally, dividing by m, we have the m-th order deformation equation where Note that according to definition (16), the right-hand side of (15) depends only on → ϑ m−1 (t). Thus, we easily gain the series ϑ 1 (t), ϑ 2 (t), . . . by solving the linear higher-order deformation Equation (15) using the well-known symbolic computation software such as Maple, Matlab, or Mathematica. Prediction and controlling the infection was studied in detail not only in [22] but also in other papers, for example [4,[23][24][25][26][27][28][29][30][31][32][33][34][35][36]. We discuss in the next section the homotopy analysis method.

Homotopy Analysis Method
For this section, we solved the interval-based uncertain model (5) with the property that where α 1 is a constant of integration. The inverse operator A −1 is given by Let the nonlinear operator be defined as The proper selection of the auxiliary parameter and function during the implementation of the HAM method can yield uniformly valid and accurate solutions [19].
Therefore, we have the m-th order deformation equation where The solution of the m-th order deformation Equation (22) for m > 1 and using h 1 = −1 and H(t) = 1 is given by Following earlier steps, we obtain and where m ≥ 1 in both last equations.

Numerical Results and Discussion
In this section, we present the results of the homotopy analysis method for solving an interval-based uncertain model. The solutions of the interval-based uncertain model with interval β = [0.01, 0.03] and constant value γ = 0.01 in Table A1, and with interval γ = [0.01, 0.015] and constant value β = 0.01 in Table A2. Tables A3 and A4 present the minimum, maximum, and midpoints of the susceptible, infected, and recovered human population with intervals of β and γ. The results of the HAM show strong agreement with the approximation technique. In Table A3, we present the result obtained by the Runge-Kutta of fourth order method for the susceptible, infected, and recovered humans. Then, we observed that the results are in good agreement with the homotopy analysis method (HAM) in Table A4.
In Table A1, we present the result of the susceptible, infected, and recovered humans, where β is considered an interval and γ is given as a constant. In Table A2, β is considered a constant and γ is given as an interval. It is observed from Table A1 that as time increases, the lower bound (minimum) and the upper bound (maximum) are decreasing for the susceptible human population. It is also detected that the lower bound (minimum) and the upper bound (maximum) of both the infected and recovered human populations increase with time.
In Table A2, it is observed that the same situation seems to be occurring in both the lower bound (minimum) and the upper bound (maximum) for the susceptible humans. It is also noticed that the lower bound (minimum) and the upper bound (maximum) of both the infected and recovered human populations increase with time.
It is seen from Tables A3 and A4 that the lower bound (minimum) and the upper bound of the susceptible population is decreasing with time, as seen from Tables A1 and A2 Tables A1 and A2, the interval [β = 0.02, γ = 0.01] represents the center for β and constant value γ. In the next section we discuss the stochastic version of the model.

Stochastic Version of the Model
In this part, we denote the complete probability space with a filtration {F t } t≥0 with (Ω, F, Q) and it satisfies the condition that it is increasing and continuous while F 0 have every Q-empty sets. We introduce randomness into Equation (5) and assume that the white noise depends on the size of the corresponding populations where we applied the corresponding pattern f i (S(t).I(t), R(t))dW(t), such that f i represents the intensity of the random perturbation i ∈ [1,3] and W(t) t≥0 is a single dimensional Brownian motion that is defined on a complete probability space Ω, F, {F t } t≥0 , Q . Then, the stochastic model of the SIR system (5) is described by the stochastic differential equations (SDEs): Let X(t) = (S(t), I(t), R(t)). Then, we can rewrite Equation (5) in the pattern of a single dimensional SDE of the form and the function G : R 2 + × R 2 + → R 2 + is given by In the next section, we discuss the graphical illustration of our results. Figure 1 shows the plot of the maximum, midpoint, and the minimum of the susceptible human intervals β = [0.01, 0.03] and γ = [0.01, 0.015]. It reveals that as the maximum value is decreasing, the midpoint is also decreasing, as is the minimum point. It is clearly seen from the plot that the uncertainty lies between the upper and lower bounds. Figure 2 shows the plot of the maximum, midpoint, and the minimum of the infected human intervals β = [0.01, 0.03] and γ = [0.01, 0.015]. It reveals that as the maximum value is increasing, the midpoint is also increasing, as is the minimum point. It is clearly seen from the plot that the uncertainty lies between the upper and lower bounds. Figure 3 shows the plot of the maximum, midpoint, and the minimum of the recovered human intervals β = [0.01, 0.03] and γ = [0.01, 0.015]. It reveals that as the maximum value is increasing, the midpoint is also increasing, as is the minimum point. It is clearly seen from the plot that the uncertainty lies between the upper and lower bounds. We discuss in Section 10 the numerical solutions of the stochastic differential equation model.

Numerical Solution of the SDE Model
In this section, we present the simulation of the SDE model (27)

Numerical Solution of the SDE Model
In this section, we present the simulation of the SDE model (27) Figure 1 shows the simulations of the dynamic behaviors of the susceptible and the infected populations under the intervals β = 0.01 and γ = 0.005. It was observed that the stochastic simulations of the susceptible and the infected populations were lower than their deterministic simulations. Figure 2 shows the simulations of the dynamic behaviors of the recovered population under the intervals β = 0.01 and γ = 0.005. It was observed that the stochastic simulations of the recovered population were higher than their deterministic simulations. Figure 3 shows the simulations of the dynamic behaviors of the susceptible and the infected populations under the intervals β = 0.03 and γ = 0.015. It was observed that the stochastic simulations of the susceptible and the infected populations were lower than their deterministic simulations. Figure 4 shows the simulations of the dynamic behaviors of the recovered populations under the intervals β = 0.03 and γ = 0.015. It was observed that the stochastic simulations of the recovered population were lower than their deterministic simulations.   Figure 2 shows the simulations of the dynamic behavior of the recovered population under the intervals = 0.01 and = 0.005. It was observed that the stochastic simulations of the recovered population were higher than their deter ministic simulations. Figure 3 shows the simulations of the dynamic behaviors of the sus ceptible and the infected populations under the intervals = 0.03 and = 0.015. It wa observed that the stochastic simulations of the susceptible and the infected population were lower than their deterministic simulations. Figure 4 shows the simulations of the dynamic behaviors of the recovered populations under the intervals = 0.03 and = 0.015. It was observed that the stochastic simulations of the recovered population were lower than their deterministic simulations.

Discussion and Conclusions
In this work, we have studied the interval-based uncertain model of a three-compart ment mathematical model rigorously. The homotopy analysis approach has been em ployed to solve the system of nonlinear equations of SIR interval uncertainty, in particular The results obtained were compared with the known solution and are found to be in good agreement. Hence, it was established here that the homotopy analysis method has greate advantages over other analytical methods in many different ways. The HAM is a serie expansion method that is directly dependent on small or large physical parameters, and hence, it is applicable for not only weakly but also strongly nonlinear problems. It also allows for the strong convergence of the solution over larger spatial and parameter do mains. It also gives excellent flexibility in the expression of the solution and how the so lution is explicitly obtained. It provides a simple way to ensure the convergence of the solution series. Comparing the stochastic and deterministic versions of the model, we saw that the population of the susceptible, infected, and recovered populations fell between the intervals obtained in the interval-based model. These suggest that the interval-based

Discussion and Conclusions
In this work, we have studied the interval-based uncertain model of a three-compartment mathematical model rigorously. The homotopy analysis approach has been employed to solve the system of nonlinear equations of SIR interval uncertainty, in particular. The results obtained were compared with the known solution and are found to be in good agreement. Hence, it was established here that the homotopy analysis method has greater advantages over other analytical methods in many different ways. The HAM is a series expansion method that is directly dependent on small or large physical parameters, and hence, it is applicable for not only weakly but also strongly nonlinear problems. It also allows for the strong convergence of the solution over larger spatial and parameter domains. It also gives excellent flexibility in the expression of the solution and how the solution is explicitly obtained. It provides a simple way to ensure the convergence of the solution series. Comparing the stochastic and deterministic versions of the model, we saw that the population of the susceptible, infected, and recovered populations fell between the intervals obtained in the interval-based model. These suggest that the interval-based model give a very good range for the general SIR epidemic model. In the future, we plan to use fuzzy differential equations to capture the dynamics, and we also plan to look into a more practical problem that may be grounded with epidemiological data.

Data Availability Statement:
The data used to support the findings of this study are included within the article.

Acknowledgments:
The authors thank the ICT Department of the Federal University Oye Ekiti, Ekiti State, Nigeria for the provision of space to complete this work and for their numerous support, as well as the Department of Science and Technology, Government of India and the Ministry of Defense of the Czech Republic for the support via the VAROPS grant.

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

Appendix A
The solutions obtained by the homotopy analysis method and the Runge-Kutta of the fourth order method for various intervals β = β, β and γ = γ, γ and for various constant values of β and γ are stated in Tables A1-A4