ScholarWorks @ UTRGV ScholarWorks @ UTRGV Integer Versus Fractional Order SEIR Deterministic and Stochastic Integer Versus Fractional Order SEIR Deterministic and Stochastic Models of Measles Models of Measles

: In this paper, we compare the performance between systems of ordinary and (Caputo) fractional differential equations depicting the susceptible-exposed-infectious-recovered (SEIR) models of diseases. In order to understand the origins of both approaches as mean-ﬁeld approximations of integer and fractional stochastic processes, we introduce the fractional differential equations (FDEs) as approximations of some type of fractional nonlinear birth and death processes. Then, we examine validity of the two approaches against empirical courses of epidemics; we ﬁt both of them to case counts of three measles epidemics that occurred during the pre-vaccination era in three different locations. While ordinary differential equations (ODEs) are commonly used to model epidemics, FDEs are more ﬂexible in ﬁtting empirical data and theoretically offer improved model predictions. The question arises whether, in practice, the beneﬁts of using FDEs over ODEs outweigh the added computational complexities. While important differences in transient dynamics were observed, the FDE only outperformed the ODE in one of out three data sets. In general, FDE modeling approaches may be worth it in situations with large reﬁned data sets and good numerical algorithms.

Fractional differential equations are usually used to involve the memory of the process in the dynamics of the systems. There is more than one type of fractional order derivative, most notably, Caputo, Grünwald-Letnikov, and Riemann-Liouville [25]. Here, we study the Caputo fractional order derivative. Integer order derivatives of ordinary differential equations are special cases of fractional order derivatives. It was noted in more than one paper, e.g., Reference [26], that FDEs give a better depiction of the courses of epidemics and natural phenomena than ODEs. Few researchers have also fitted their FDE models to data, e.g., Reference [26,27]. This motivated us to compare systems of ODEs and FDEs by fitting them to some actual epidemic data.
Measles is a marker disease for virological, epidemiological, clinical, statistical, geographical, mathematical, and humanitarian reasons [28] (pp. [16][17][18][19][20][21]. Mathematical modeling of measles epidemics dates back as far as 1888 by D'Enko and then by Hamer; see Haggett [28] (p. 19). Regularity and a large number of cases of measles epidemics with major peaks in the pre-vaccination era (before 1964) support the choice of testing models on measles data. Many other researchers formulated measles models and fit them to data, as in Bjørnstad et al. [29], where a time scale of two weeks is recommended fitting the number of cases, and in Yingcun Xia et al. [30], where a model is used to examine a spatial network. In this paper, we choose to use data of measles infections in the USA and UK in two decades of the pre-vaccination era  to compare the goodness of fit of ODEs and FDEs to those epidemics.
While ODEs are well-established as deterministic models of the spread of diseases [31,32], FDE models are sometimes used. However, often these approaches lack mathematical basis or physical interpretation except for exchanging integer differentiation with fractional ones [26,33]. Angstmann et al. [34] and Sardar et al. [35] provided a valid variation by considering the memory of the non-Markovian infection process. The result is a mixed system of integer and fractional derivatives of the Riemann-Liouville type. Saeedian et al. [36] showed how another memory functional of the process can lead to replacing the integer derivatives with Caputo fractional derivatives. In this paper, we show how Caputo fractional differential equations follow naturally from fractional stochastic processes like those introduced in literature [37][38][39][40][41][42][43][44][45][46]. We then compare transient and long term dynamics between the FDE and ODE models while fitting them to three different data sets. The Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) are used to compare between the fits of the two models to three data sets. For completeness, we will cover all the required background and the relevant definitions in Section 2. That includes a synopsis of Caputo's fractional calculus and fractional stochastic SEIR processes. Section 2 will also include the derivation of the fractional order differential equation depicting the SEIR model from the fractional stochastic SEIR process. It will be followed by the stability analysis of the equilibria of the system of fractional differential equations, which will be then fitted to measles data and simulated.

Methods
In this section, we provide a background for fractional differentiation and a fractional birth and death process. We also introduce the integer and fractional differential equations for the SEIR model and analyze the stability of the FDEs' equilibria.

Fractional Calculus
Let D n be the Leibniz integer-order differential operator given by and let J n be an integration operator of integer order given by where n ∈ Z + . Let us use D = D 1 for the first derivative. We will use ∂ α x F := ∂ α F ∂x α and use ∂ x F := ∂F ∂x .

Fractional Stochastic Process
Fix 0 < α ≤ 1. Following Earn et al. [49], we consider a compartmental susceptible-exposedinfected-recovered (SEIR) model to depict the measles transmission dynamics in a closed population.
be the number of susceptible, exposed, infected, and recovered individuals, respectively, such that X = N, the population size. Figure 1 shows how the disease is progressing from one sub-population to another. Here, µ is the recruitment and per capita death rate, β is the transmission rate, δ is the rate at which exposed individuals become infectious, and σ is the recovery rate.
A stochastic SEIR model can be depicted using a continuous time Markov chain (CTMC) like the birth and death process with non-linear rates of transition as those given in Table 1; see Allen [50] (p. 22) and [51] (p. 321). Bartlett, M.S. [9] and Greenwood and Gordillo [31] introduced (integer) stochastic SIR model using CTMC with rates similar to those in the first six rows in Table 1 to show a deterministic SIR model of the ODE type depicting the approximate dynamics of the means of the processes. Here, we introduce a fractional SEIR model using a CTMC of fractional birth and death process on triplets (i, j, k) with rates provided by Table 1. Table 1. Transitions and their rates for a birth and death process depicting a stochastic SEIR model.

Transition
Rate An α-fractional SEIR stochastic process {(X for i, j, k = 0, 1, . . ., such that 0 ≤ i + j + k ≤ N and P((X (α) 3 (0)) = (i 0 , j 0 , k 0 )) = 1, has a fractional forward Kolmogorov equation of the stochastic SEIR model similar to Equation (A1) and is given by with p (α) (i,j,k) (t) = 0 if either i, j, or k are negative or i + j + k > N (see Di Crescenzo et al. [45]). The classical forward Kolmogorov equation of the stochastic SEIR model follows when α = 1 with state probabilities p (1) (i,j,k) (t), [51] (p. 321). Equation (10) can be used to find the probability generating function ) of the state probabilities, as the solution of the Cauchy problem for t > 0 and Note that the integer or classical stochastic SEIR process is (X (1) 3 (t)) which is simply the case when α = 1. Interestingly, the fractional stochastic SEIR process is a random-time subordination of the integer stochastic SEIR model, as established for other fractional processes like the fractional Poisson process [37,45,52], and the fractional birth and/or death processes [39,40,42,43,53]. In Mandelbrot and Taylor [54], the stochastic time process T 2α is called the operational time, and t is the physical time.
3 (t)) has the same distribution as the random-time subordinated integer stochastic SEIR process The proof is provided in Appendix A.

Measles Model via Fractional Differential Equations (FDE)
The means of the three discrete-marginal processes X , where N is the total population size and E(x) is the expected value of x. Thus, using Equation (11) and approximating 3 (t)), we reach the fractional order version of the system of equations that was used by Bartlett, M.S. [9] to model measles: where S (α) , E (α) , and I (α) be the proportion of susceptible, exposed, and infected individuals, respectively. With proportion of recovered individuals given by R (α) = 1 − (S (α) + E (α) + I (α) ), we reach the fractional α order SEIR model with 0 < α ≤ 1. The non-negative parameters β, µ, δ, and σ-denoting them by θ, for brevity-have dimensions given by 1 time α . By construction of the FDE model as a mean field approximation of the α-fractional stochastic SEIR process which, in its turn, is a subordination of an integer stochastic SEIR process by Theorem 1, those parameters could be interpreted as the rates measured by an independent observer of the process or calculated based on a cosmic time flow [47]. We replace those parameters with a power α of new parameters; that is, θ α * in place of θ so the parameters θ * will have the dimension of 1 time , and the system becomes the following form: By replacing θ by θ α * , the system of integer order differential equations with its epidemiological parameters follows directly from the system of fractional order differential equations when α = 1.

Measles Model via Ordinary Differential Equations (ODE)
The following system of differential equations represents the ordinary differential equation representation of the SEIR model and is the FDE model when α = 1 in Equation (14).
where µ, β, δ, and σ are the model parameters described above. They all have dimensions given by 1 time . The last Equation in (15) is redundant since R = 1 − (S + E + I).

Measles Model via α-Dependent Ordinary Differential Equations
We are interested in comparing the FDE versus ODE modeling approaches. It is important to note that the basic ODE case considers α = 1; however, in the FDE case, α appears in the derivative, as well as the parameter, values. In order to better compare these two approaches, here, we develop an ODE analogue to the FDE that incorporates α in the parameter values. We call this new system the α-dependent ODE. By dropping the α order derivative from the left side and α power from S (α) , E (α) , and I (α) of Equation (14), our α-dependent ODE takes the following form: Here, the above systems may differ from the classical ODE (Equation (15)) if α = 1. In this case, the α value used in the above α-dependent ODE are obtained from data fitting procedures using the FDE model (Equation (14)).

Model Analysis
Analysis of the ODE is almost the same as of the FDE, so we include the FDE one here. We start by proving the positive invariance of the region of solutions of the FDE model. Henceforth, we drop the α from S (α) , E (α) , and I (α) , for brevity.
The following two lemmas of asymptotic behavior of FDEs are given here and their proof in Appendix A for completeness.
is a positive invariant set for the FDE model in Equation (14).
Thus, in general, the stability of the ordinary differential equations model implies stability of its fractional counter model of order α. But, here, they are equivalent due to the following lemma, in which the solution could be found in Appendix A.

Lemma 2.
The Disease-Free Equilibrium DFE is locally asymptotically stable if R 0 < 1. The endemic equilibrium EE is is locally asymptotically stable if R 0 > 1. Therefore, they have the same asymptotic behavior. Yet, the transient behavior differs, as will be seen by simulations below.
Moreover, a very important difference is their oscillation behavior is not similar. Let λ and u for = 1, 2, . . . , n be the eigenvalues and their respective eigenvectors of an n × n matrix A. The general solution of initial value problem consisting of a system of n linear fractional differential equations D α * x(t) = Ax(t), such that x(0) = x 0 , can be found to be for certain constants c ∈ C for = 1, 2, . . . , n, such that ∑ n =1 c u = x 0 , [25] (Theorem 7.13). In case that α = 1, we recover the known solution of the system of ODEs given by If n = 3 and A is not a symmetric matrix, then at least one of the eigenvalues is a real-valued number and the other two eigenvalues , say λ 2 and λ 3 , are conjugate complex-valued. In that situation, x(t) would oscillate with inter-peak periods, called inter-epidemic period in disease modeling, given by 2π(Im(λ 2 )) −1 [14]. If Re(λ ) < 0 for all , then the oscillations will be damped to zero. That damped oscillation is straightforward in the case of α = 1 due to the exponential damping in the superposition of the sine and cosine functions. That behavior, however, is not straightforward for 0 < α < 1.

Numerical Simulations
Since the mean of the subordinator process is E(T α (t)) = t α Γ(α + 1) , we use a method similar to that which was introduced in Demirici and Özalp [55] to find approximate solutions to initial value FDE problems. We use that method here to simulate the solution of the FDE measles SEIR model. Consider the initial value problem: for some T > 0. A solution of Equation (19) is approximated by the deterministic time subordination of y(s), the solution of the ordinary differential equation dy(s) ds = g(s, y(s)), for 0 < s ≤ t α Γ(α + 1) where for all 0 < t ≤ T, [55].
We use the subordination of the solution of ODEs to FDEs represented in Equations (20) and (21) to numerically simulate solutions of FDEs; see Algorithm 1.

Fitting FDE and ODE Models to Measles Data
We use the method of ordinary least squares (OLS) to fit the FDE model to the data by minimizing the objective function for α ∈ (0, 1], and β, µ, δ, σ ∈ (0, ∞), where {I i , i = 1, . . . , n} is the data of actual proportion of infections and {Î i , i = 1, . . . , n} is the simulated proportion of infections. The valuesÎ i approximating I(t i ) are found by solving the FDE model using Algorithm 1.
Parameter estimation was conducted using MATLAB MultiStart and fmincon functions. The MultiStart function carries out the optimization procedure using initial points within the parameters' spaces. It generates some initial points depending on a converging algorithm. The fmincon function finds a local minimum for the constrained nonlinear multivariable function. The MultiStart, together with fmincon functions, does the global optimization of a nonlinear multivariable function. The MultiStart function uses parallel processing, which drastically reduces the running time.

Results
We solve the system of FDE (Equation (14)) using Algorithm 1 and the systems of ODE (Equations (15) and (16)) using the Runge-Kutta method.
In order to study the qualitative aspects of the FDE and its numerical solution, we carried out a number of simulations. Simulations of the classical ODE (Equation (15)) and FDE (Equation (14)), Figure 2, shows that the system of fractional differential equations is very sensitive to its order of differentiation α. For smaller α, the peak number of cases of the epidemic is larger but the duration of the outbreak is shorter. The solution of the FDE model converges to the solution of the classical ODE as α → 1. To further compare the two modeling approaches, we consider the analogue ODEs derived for specific α values; see Equation (16). These comparisons are shown in Figure 3. During transient dynamics both models exhibit several peaks in the number of cases. The number of these peaks and their respective amplitudes are similar between models; however, there are differences in the timing of these peaks. The transient oscillations of the FDE model are more stretched out than its ODE analogue, and its solutions experience longer inter-epidemic times. Both models approach the same equilibria solutions.
Simulations of Equation (18), Figure 4, shows that disease models of fractional order equations lack the same oscillatory behavior exhibited by systems of ODEs with conjugate complex eigenvalues of the Jacobian matrices calculated at endemic equilibrium.   The models were fitted to three measles epidemics in the pre-vaccination era in three different cities: New York, London, and Portsmouth. Simulations of the fitted ODE and FDE models are shown in Figure 5. See Appendix B for the data and the parameter estimates, as well.  (15)). Statistical measurements on the data fittings are presented in Table A4, found in Appendix B.4. The Akaike information criterion (AIC), Bayesian information criterion (BIC), and Sum squared error (SSE) are very similar between the two models for all three data sets. However, the biggest differences between the two models are found with the Portsmouth data fittings where the FDE performs slightly better than the ODE.

Discussion
Replacing first order derivatives with Caputo fractional derivatives has been the practice for many studies using fractional order modeling of diseases. In this paper, we show how these models follow from an approximation to the dynamical system governing the means of fractional stochastic SEIR processes. Moreover, we study ordinary and fractional order systems of differential equations of SEIR models using three data sets of measles epidemics in three different cities selected from the pre-vaccination era. Two of the data sets fits resulted in α values of extremely close to one, with differences appearing in the ninth decimal digit. In these two cases, the two models are equivalent, apart from round off errors potentially from the numerical methods. We suspect these round off error are also causing the extremely slight differences in the statistical measures (AIC, BIC, SSE) of the data fittings. In the third case, the data for Portsmouth, where α = 0.88 the FDE model outperformed the ODE model; see Figure 5 and Appendix B.4.
Angstmann et al. [34] use the master equation of a continuous-time random walk to derive an FDEM involving Riemann-Liouville fractional derivatives. Power laws are postulated to model time of infectiousness and recovery. That extension from exponential times in ordinary differential equations is a different approach from the mean field approximation of a stochastic process. Saeedian et al. [36] introduced the Caputo fractional differential equations through a memory of the whole process of infection and disease recovery. In our paper, we have considered, for the first time, fractional stochastic SEIR model and have shown how the Caputo fractional differential equations follows as mean-field approximation of the process. Fractional stochastic SEIR model introduced here turns out to be a random-time subordination of a classical stochastic SEIR model. Other real-life systems are modeled using a subordination of a stochastic process. A subordinated process was introduced by Mandelbrot and Taylor [54] to model the logarithm of market prices where the original process is a Brownian motion subordinated by a stochastic time process T 2α , which is the same random time process we have found here.
Further study of the fractional stochastic SEIR model might lead to interesting dynamical behaviors. For instance, it can provide more insights into the stochastic oscillations of the disease in a more flexible way than their classical counterparts. Thus, studying the fractional stochastic SEIR model is the next step in this work.

Conclusions
In this paper, we compare two deterministic models of disease: ordinary differential equations (ODE) and fractional differential equations (FDE). We use three different data sets of measles epidemics from the pre-vaccination era. We also explain FDEs as the mean-field approximation of a fractional stochastic SEIR model. To our knowledge, this is the first time such a fractional stochastic process is introduced and connected to the fractional order differential equations.
We conclude that, depending on the specific data set, the FDE model either converges to the ODE model and fits data similarly, or fits the data better and outperforms the ODE model. It is worth asking whether the complexities introduced by FDE approach are worth pursuing, as the numerical algorithms for solving FDEs are computationally expensive. Here, the FDE approach yielded better results in only one of three data sets, and the results were still similar to those obtained by the ODE. However, other data sets may yield larger differences between these two modeling approaches and benefit more from an FDE approach. In particular, more refined and larger data sets may be more appropriate for this approach. Here, oscillations can persist longer in the FDE solutions than those of the ODE (Figure 3, α = 0.85), suggesting perhaps refined data set with seasonal variations may be good candidates for the FDE modeling approach. It is also important to note that the numerical methods for solving FDEs and data fitting procedures can be improved. Future studies, with good empirical data sets, should consider the FDE approach to modeling epidemics, as well as improving numerical algorithms. The inverse transform is defined by where C is a contour parallel to the imaginary axis and to the right of the singularities off . The Laplace transform of the Caputo fractional derivative is given by

Fractional Birth and Death Process
An α-fractional nonlinear birth and death process {N α (t) : t ≥ 0} for 0 < α ≤ 1 with state probabilities p α n (t) = P(N α (t) = n|N α (0) = 1) for n ≥ 0 is defined through the forward Kolmogorov (difference-)differential equations for n ≥ 0 [39,43,53]. The rates λ n and µ n are non-negative. The classical birth and death process follows when α = 1 with state probabilities p 1 n (t). When λ n = λ and µ n = 0 for all n, the α-fractional nonlinear birth and death process becomes the α-fractional Poisson process [37,41,46]. There, it has shown that N α (t) has the same probability distribution as N(T α (t)), where N(t) is the classical birth and death process which is independent of a random time process T α (t); that is, a birth and death process subordinated by an α-stable time process.
Note that, the absolute values of partial derivatives of G are finite; that is, |∂ (i,j,k) u,v,w G| < ∞ for any i, j, k = 0, 1, 2, . . .. That is true since |u|, |v|, |w| < 1 and the population size is finite. Thus, switching integrals with derivatives or summations below are valid.
Proof of Theorem 1. We are going to show that Laplace transform of the probability generating function of the process (X (1) 3 (T 2α (t))) is the same as Laplace transformĜ of G, that solves Equation (11). From there we will conclude that the two probability distributions are the same since the probability generating function of (X (α) 3 (t)), by construction, is also a solution to the Cauchy problem in Equation (11). From Equation (11), the Laplace transformĜ is the solution of Let H (α) (u, v, w, t) be the probability generating function of the state probabilities q (α) (i,j,k) (t) = P((X 3 (T 2α (0))) = (i 0 , j 0 , k 0 )). That means that Thus, the Laplace transform of the probability generating function H (α) is given bŷ Now, the Laplace transform of the probability generating function of the process (X (1) 3 (t)) also solves (A2) when α = 1 which is If we substitute with s = r α in Equation (A3) and multiply both sides by r α−1 we get which is the same as Equation (A2). This completes the proof.