Integer Versus Fractional Order SEIR Deterministic and Stochastic Models of Measles

Abstract

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-field 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 fit 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 flexible in fitting empirical data and theoretically offer improved model predictions. The question arises whether, in practice, the benefits 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 refined data sets and good numerical algorithms.

1. Introduction

Modeling the spread of infectious diseases before the introduction of vaccines, as well as the validation of these models, has been widely studied since the works of Reference [1,2,3,4,5,6,7,8,9,10]. See Bailey [11] and Anderson [12] for more details about the history of disease modeling. Deterministic models using ordinary differential equations (ODEs) have received great attention [12,13,14,15,16] and wide assimilation by health sciences [17]. Other deterministic models, such as difference equations, are also used to model the spread of diseases; for instance, see Fisman et al. [18]. However, fractional differential equations (FDEs) have been used in the last decade to model the course of epidemics [19,20,21,22,23,24].

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–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 (1944–1964) 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.

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

2.1. Preliminaries

Fractional Calculus

Let $D^n$ be the Leibniz integer-order differential operator given by

$$D^{n}f=\frac{d^{n}f}{d{t}^n}=f^{\left(n\right)},$$

and let $J^n$ be an integration operator of integer order given by

$$J^{n}f\left(t\right)=\frac{1}{n!}{\int }_0^t{(t-\tau )}^{n-1}f\left(\tau \right)d\tau ,$$

(1)

where $n\in {\mathbb{Z}}^{+}$. Let us use $D=D^1$ for the first derivative. We will use ${\partial }_x^{\mathsf{\alpha }}F:={\displaystyle \frac{{\partial }^{\mathsf{\alpha }}F}{\partial x^{\mathsf{\alpha }}}}$ and use ${\partial }_{x}F:={\displaystyle \frac{\partial F}{\partial x}}$.

For fraction-order integrals, we use

$$J^{n-\mathsf{\alpha }}f\left(t\right)=\frac{1}{\mathrm{\Gamma }(n-\mathsf{\alpha })}{\int }_0^t{(t-\tau )}^{n-\mathsf{\alpha }-1}f\left(\tau \right)d\tau ,$$

(2)

where $n-1<\mathsf{\alpha }\le n$. Now, define the Caputo fractional differential operator $D_{\ast }^{\mathsf{\alpha }}$ to be

$$D_{\ast }^{\mathsf{\alpha }}f\left(t\right)=J^{n-\mathsf{\alpha }}{D}^{n}f\left(t\right),$$

where $n-1<\mathsf{\alpha }\le n$, for $n\in \mathbb{N}$. It is also known that

$$\begin{array}{cc}\hfill \underset{\mathsf{\alpha }\to n}{lim}{D}_{\ast }^{\mathsf{\alpha }}f\left(t\right)& =f^{\left(n\right)}\left(t\right),\hfill \\ \hfill \underset{\mathsf{\alpha }\to n-1}{lim}{D}_{\ast }^{\mathsf{\alpha }}f\left(t\right)& =f^{(n-1)}\left(t\right)-f^{(n-1)}\left(0\right)\hfill \end{array}$$

(3)

for any $n\in \mathbb{N}$. We will consider $n=1$ in this work; that is $0<\mathsf{\alpha }\le 1$. In that case,

$$J^{1-\mathsf{\alpha }}f\left(t\right)={\int }_0^{t}f\left(\tau \right)d{g}_t\left(\tau \right),$$

(4)

where $g_t\left(\tau \right)=\frac{1}{\mathrm{\Gamma }(2-\mathsf{\alpha })}\left(t^{1-\mathsf{\alpha }}-{(t-\tau )}^{1-\mathsf{\alpha }}\right)$. That is, for each t, the integral $J^{1-\mathsf{\alpha }}f\left(t\right)$ is an area under $f\left(\tau \right)$, while above the curve of $g_t\left(\tau \right)$ that works as a deformed or slowed time-scale, as illustrated by Podlubny [47].

The generalized mean-value theorem for the Caputo fractional derivative is given as

$$f\left(x\right)=f\left(a\right)+\frac{1}{\mathrm{\Gamma }\left(\mathsf{\alpha }\right)}{D}_{\ast }^{\mathsf{\alpha }}f\left(c\right){(x-a)}^{\mathsf{\alpha }}\mathrm{for}\mathrm{some}a\le c\le x$$

and for all $x\in (a,b]$ whenever $f,D_{\ast }^{\mathsf{\alpha }}f\in C\left([a,b]\right)$; see, e.g., Özalp and Demirci [48].

The Mittag-Leffler is a function that generalizes the exponential function. That function can be written as follows:

$$E_{\mathsf{\alpha }}\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathrm{\Gamma }(\mathsf{\alpha }k+1)},\mathsf{\alpha }\in {\mathbb{R}}^{+},z\in \mathbb{C},$$

(5)

or, more generally using two parameters,

$$E_{\mathsf{\alpha },\mathsf{\beta }}\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathrm{\Gamma }(\mathsf{\alpha }k+\mathsf{\beta })},\mathsf{\alpha },\mathsf{\beta }\in {\mathbb{R}}^{+},z\in \mathbb{C}.$$

(6)

The general Mittag-Leffler has the following important property for any $\mathsf{\alpha },\mathsf{\beta }>0$:

$$E_{\mathsf{\alpha },\mathsf{\beta }}\left(z\right)=z{E}_{\mathsf{\alpha },\mathsf{\alpha }+\mathsf{\beta }}\left(z\right)+\frac{1}{\mathrm{\Gamma }\left(\mathsf{\beta }\right)}.$$

(7)

Two important differential properties of the Mittag-Leffler function are that

$$D_{\ast }^{\mathsf{\alpha }}{e}^{\mathsf{\lambda }t}=t^{-\mathsf{\alpha }}{E}_{1,1-\mathsf{\alpha }}\left(\mathsf{\lambda }t\right)$$

(8)

and

$$D_{\ast }^{\mathsf{\alpha }}{E}_{\mathsf{\alpha },1}\left(\mathsf{\lambda }{t}^{\mathsf{\alpha }}\right)=\mathsf{\lambda }{E}_{\mathsf{\alpha },1}\left(\mathsf{\lambda }{t}^{\mathsf{\alpha }}\right)$$

(9)

for any $\mathsf{\lambda }>0$.

2.2. Fractional Stochastic Process

Fix $0<\mathsf{\alpha }\le 1$. Following Earn et al. [49], we consider a compartmental susceptible-exposed- infected-recovered (SEIR) model to depict the measles transmission dynamics in a closed population. Let $X_1^{\left(\mathsf{\alpha }\right)}$, $X_2^{\left(\mathsf{\alpha }\right)}$, $X_3^{\left(\mathsf{\alpha }\right)}$, and $X_4^{\left(\mathsf{\alpha }\right)}$ be the number of susceptible, exposed, infected, and recovered individuals, respectively, such that $X_1^{\left(\mathsf{\alpha }\right)}+X_2^{\left(\mathsf{\alpha }\right)}+X_3^{\left(\mathsf{\alpha }\right)}+X_4^{\left(\mathsf{\alpha }\right)}=N$, the population size. Figure 1 shows how the disease is progressing from one sub-population to another.

Here, $\mathsf{\mu }$ is the recruitment and per capita death rate, $\mathsf{\beta }$ is the transmission rate, $\mathsf{\delta }$ is the rate at which exposed individuals become infectious, and $\mathsf{\sigma }$ 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. Transitions and their rates for a birth and death process depicting a stochastic SEIR model.

Transition	Rate
$X_1^{\left(\mathsf{\alpha }\right)}\to X_1^{\left(\mathsf{\alpha }\right)}+1$	$\mathsf{\mu }N$
$X_1^{\left(\mathsf{\alpha }\right)}\to X_1^{\left(\mathsf{\alpha }\right)}-1$	$\mathsf{\beta }{X}_1^{\left(\mathsf{\alpha }\right)}{\displaystyle \frac{X_3^{\left(\mathsf{\alpha }\right)}}{N}}+\mathsf{\mu }{X}_1^{\left(\mathsf{\alpha }\right)}$
$X_2^{\left(\mathsf{\alpha }\right)}\to X_2^{\left(\mathsf{\alpha }\right)}+1$	$\mathsf{\beta }{X}_1^{\left(\mathsf{\alpha }\right)}{\displaystyle \frac{X_3^{\left(\mathsf{\alpha }\right)}}{N}}$
$X_2^{\left(\mathsf{\alpha }\right)}\to X_2^{\left(\mathsf{\alpha }\right)}-1$	$(\mathsf{\mu }+\mathsf{\delta })X_2^{\left(\mathsf{\alpha }\right)}$
$X_3^{\left(\mathsf{\alpha }\right)}\to X_3^{\left(\mathsf{\alpha }\right)}+1$	$\mathsf{\delta }{X}_2^{\left(\mathsf{\alpha }\right)}$
$X_3^{\left(\mathsf{\alpha }\right)}\to X_3^{\left(\mathsf{\alpha }\right)}-1$	$(\mathsf{\mu }+\mathsf{\sigma })X_3^{\left(\mathsf{\alpha }\right)}$
$X_4^{\left(\mathsf{\alpha }\right)}\to X_4^{\left(\mathsf{\alpha }\right)}+1$	$\mathsf{\sigma }{X}_3^{\left(\mathsf{\alpha }\right)}$
$X_4^{\left(\mathsf{\alpha }\right)}\to X_4^{\left(\mathsf{\alpha }\right)}-1$	$\mathsf{\mu }{X}_4^{\left(\mathsf{\alpha }\right)}$

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

An $\mathsf{\alpha }$-fractional SEIR stochastic process $\{(X_1^{\left(\mathsf{\alpha }\right)}\left(t\right),X_2^{\left(\mathsf{\alpha }\right)}\left(t\right),X_3^{\left(\mathsf{\alpha }\right)}\left(t\right)):t\ge 0\}$ for $0<\mathsf{\alpha }\le 1$ with state probabilities

$$p_{(i,j,k)}^{\left(\mathsf{\alpha }\right)}\left(t\right)=P((X_1^{\left(\mathsf{\alpha }\right)}\left(t\right),X_2^{\left(\mathsf{\alpha }\right)}\left(t\right),X_3^{\left(\mathsf{\alpha }\right)}\left(t\right))=(i,j,k)|(X_1^{\left(\mathsf{\alpha }\right)}\left(0\right),X_2^{\left(\mathsf{\alpha }\right)}\left(0\right),X_3^{\left(\mathsf{\alpha }\right)}\left(0\right))=(i_0,j_0,k_0))$$

for $i,j,k=0,1,\dots$, such that $0\le i+j+k\le N$ and $P((X_1^{\left(\mathsf{\alpha }\right)}\left(0\right),X_2^{\left(\mathsf{\alpha }\right)}\left(0\right),X_3^{\left(\mathsf{\alpha }\right)}\left(0\right))=(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

$$\begin{array}{ll}{D}_{\ast }^{\mathsf{\alpha }}{p}_{(i,j,k)}^{\left(\mathsf{\alpha }\right)}\left(t\right)& =\mathsf{\mu }N{p}_{(i-1,j,k)}^{\left(\mathsf{\alpha }\right)}\left(t\right)+\mathsf{\beta }(i+1)\frac{k}{N}{p}_{(i+1,j-1,k)}^{\left(\mathsf{\alpha }\right)}\left(t\right)+\mathsf{\mu }(i+1)p_{(i+1,j,k)}^{\left(\mathsf{\alpha }\right)}\left(t\right)\\ & +\mathsf{\delta }(j+1)p_{(i,j+1,k-1)}^{\left(\mathsf{\alpha }\right)}\left(t\right)+\mathsf{\mu }(j+1)p_{(i,j+1,k)}^{\left(\mathsf{\alpha }\right)}\left(t\right)+(\mathsf{\sigma }+\mathsf{\mu })(k+1)p_{(i,j,k+1)}^{\left(\mathsf{\alpha }\right)}\left(t\right)\\ & -(\mathsf{\mu }N+\mathsf{\beta }i\frac{k}{N}+\mathsf{\mu }i+(\mathsf{\delta }+\mathsf{\mu })j+(\mathsf{\sigma }+\mathsf{\mu })k)p_{(i,j,k)}^{\left(\mathsf{\alpha }\right)}\left(t\right)\end{array},$$

(10)

with $p_{(i,j,k)}^{\left(\mathsf{\alpha }\right)}\left(t\right)=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 $\mathsf{\alpha }=1$ with state probabilities $p_{(i,j,k)}^{\left(1\right)}\left(t\right)$,51] (p. 321). Equation (10) can be used to find the probability generating function $G^{\left(\mathsf{\alpha }\right)}(u,v,w,t)=E\left(u^{X_1^{\left(\mathsf{\alpha }\right)}\left(t\right)}{v}^{X_2^{\left(\mathsf{\alpha }\right)}\left(t\right)}{w}^{X_3^{\left(\mathsf{\alpha }\right)}\left(t\right)}\right)$ of the state probabilities, as the solution of the Cauchy problem

$$\begin{array}{ll}{D}_{\ast }^{\mathsf{\alpha }}{G}^{\left(\mathsf{\alpha }\right)}& =\mathsf{\mu }N(u-1)G^{\left(\mathsf{\alpha }\right)}+\mathsf{\mu }(1-u){\partial }_u{G}^{\left(\mathsf{\alpha }\right)}+(\mathsf{\delta }w+\mathsf{\mu }-(\mathsf{\delta }+\mathsf{\mu })v){\partial }_v{G}^{\left(\mathsf{\alpha }\right)}\\ & +(\mathsf{\sigma }+\mathsf{\mu })(1-w){\partial }_w{G}^{\left(\mathsf{\alpha }\right)}+\mathsf{\beta }\frac{w}{N}(v-u){\partial }_{uw}{G}^{\left(\mathsf{\alpha }\right)}\end{array}$$

(11)

for $t>0$ and $G^{\left(\mathsf{\alpha }\right)}(u,v,w,0)=u^{i_0}{v}^{j_0}{w}^{k_0}$, for $-1<u,v,w<1$.

Note that the integer or classical stochastic SEIR process is $(X_1^{\left(1\right)}\left(t\right),X_2^{\left(1\right)}\left(t\right),X_3^{\left(1\right)}\left(t\right))$ which is simply the case when $\mathsf{\alpha }=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 ${\mathcal{T}}_{2\mathsf{\alpha }}$ is called the operational time, and t is the physical time.

Theorem 1.

The fractional stochastic SEIR process $(X_1^{\left(\mathsf{\alpha }\right)}\left(t\right),X_2^{\left(\mathsf{\alpha }\right)}\left(t\right),X_3^{\left(\mathsf{\alpha }\right)}\left(t\right))$ has the same distribution as the random-time subordinated integer stochastic SEIR process

$$(X_1^{\left(1\right)}\left({\mathcal{T}}_{2\mathsf{\alpha }}\left(t\right)\right),X_2^{\left(1\right)}\left({\mathcal{T}}_{2\mathsf{\alpha }}\left(t\right)\right),X_3^{\left(1\right)}\left({\mathcal{T}}_{2\mathsf{\alpha }}\left(t\right)\right))$$

for $t>0$ and $0<\mathsf{\alpha }\le 1$.

The proof is provided in Appendix A.

2.3. Measles Model via Fractional Differential Equations (FDE)

The means of the three discrete-marginal processes $X_1^{\left(\mathsf{\alpha }\right)}\left(t\right)$, $X_2^{\left(\mathsf{\alpha }\right)}\left(t\right)$, and $X_3^{\left(\mathsf{\alpha }\right)}\left(t\right)$ can be found using ${\partial }_u{G}^{\left(\mathsf{\alpha }\right)}(1,1,1,t)$, ${\partial }_v{G}^{\left(\mathsf{\alpha }\right)}(1,1,1,t)$, and ${\partial }_w{G}^{\left(\mathsf{\alpha }\right)}(1,1,1,t)$, respectively. Let $S^{\left(\mathsf{\alpha }\right)}\left(t\right):=\frac{1}{N}\mathbb{E}(X_1^{\left(\mathsf{\alpha }\right)}\left(t\right))$, $E^{\left(\mathsf{\alpha }\right)}\left(t\right):=\frac{1}{N}\mathbb{E}(X_2^{\left(\mathsf{\alpha }\right)}\left(t\right))$, and $I^{\left(\mathsf{\alpha }\right)}\left(t\right):=\frac{1}{N}\mathbb{E}(X_3^{\left(\mathsf{\alpha }\right)}\left(t\right))$, where N is the total population size and $\mathbb{E}\left(x\right)$ is the expected value of x. Thus, using Equation (11) and approximating $\mathbb{E}(X_1^{\left(\mathsf{\alpha }\right)}\left(t\right)X_3^{\left(\mathsf{\alpha }\right)}\left(t\right))$ by $\mathbb{E}(X_1^{\left(\mathsf{\alpha }\right)}\left(t\right))E(X_3^{\left(\mathsf{\alpha }\right)}\left(t\right))$, we reach the fractional order version of the system of equations that was used by Bartlett, M.S. [9] to model measles:

$$\begin{array}{cc}\hfill D_{\ast }^{\mathsf{\alpha }}{S}^{\left(\mathsf{\alpha }\right)}& =\mathsf{\mu }-\mathsf{\beta }{S}^{\left(\mathsf{\alpha }\right)}{I}^{\left(\mathsf{\alpha }\right)}-\mathsf{\mu }{S}^{\left(\mathsf{\alpha }\right)}\hfill \\ \hfill D_{\ast }^{\mathsf{\alpha }}{E}^{\left(\mathsf{\alpha }\right)}& =\mathsf{\beta }{S}^{\left(\mathsf{\alpha }\right)}{I}^{\left(\mathsf{\alpha }\right)}-(\mathsf{\mu }+\mathsf{\delta })E^{\left(\mathsf{\alpha }\right)}\hfill \\ \hfill D_{\ast }^{\mathsf{\alpha }}{I}^{\left(\mathsf{\alpha }\right)}& =\mathsf{\delta }{E}^{\left(\mathsf{\alpha }\right)}-(\mathsf{\mu }+\mathsf{\sigma })I^{\left(\mathsf{\alpha }\right)}\hfill \end{array}$$

(12)

where $S^{\left(\mathsf{\alpha }\right)}$, $E^{\left(\mathsf{\alpha }\right)}$, and $I^{\left(\mathsf{\alpha }\right)}$ be the proportion of susceptible, exposed, and infected individuals, respectively. With proportion of recovered individuals given by $R^{\left(\mathsf{\alpha }\right)}=1-(S^{\left(\mathsf{\alpha }\right)}+E^{\left(\mathsf{\alpha }\right)}+I^{\left(\mathsf{\alpha }\right)})$, we reach the fractional $\mathsf{\alpha }$ order SEIR model

$$\begin{array}{cc}\hfill D_{\ast }^{\mathsf{\alpha }}{S}^{\left(\mathsf{\alpha }\right)}& =\mathsf{\mu }-\mathsf{\beta }{S}^{\left(\mathsf{\alpha }\right)}{I}^{\left(\mathsf{\alpha }\right)}-\mathsf{\mu }{S}^{\left(\mathsf{\alpha }\right)}\hfill \\ \hfill D_{\ast }^{\mathsf{\alpha }}{E}^{\left(\mathsf{\alpha }\right)}& =\mathsf{\beta }{S}^{\left(\mathsf{\alpha }\right)}{I}^{\left(\mathsf{\alpha }\right)}-(\mathsf{\mu }+\mathsf{\delta })E^{\left(\mathsf{\alpha }\right)}\hfill \\ \hfill D_{\ast }^{\mathsf{\alpha }}{I}^{\left(\mathsf{\alpha }\right)}& =\mathsf{\delta }{E}^{\left(\mathsf{\alpha }\right)}-(\mathsf{\mu }+\mathsf{\sigma })I^{\left(\mathsf{\alpha }\right)}\hfill \\ \hfill D_{\ast }^{\mathsf{\alpha }}{R}^{\left(\mathsf{\alpha }\right)}& =\mathsf{\sigma }{I}^{\left(\mathsf{\alpha }\right)}-\mathsf{\mu }{R}^{\left(\mathsf{\alpha }\right)}\hfill \end{array}$$

(13)

with $0<\mathsf{\alpha }\le 1$. The non-negative parameters $\mathsf{\beta }$, $\mathsf{\mu }$, $\mathsf{\delta }$, and $\mathsf{\sigma }$—denoting them by $\theta$, for brevity—have dimensions given by $\frac{1}{{\mathrm{time}}^{\mathsf{\alpha }}}$. By construction of the FDE model as a mean field approximation of the $\mathsf{\alpha }$-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 $\mathsf{\alpha }$ of new parameters; that is, ${\theta }_{\ast }^{\mathsf{\alpha }}$ in place of $\theta$ so the parameters ${\theta }_{\ast }$ will have the dimension of $\frac{1}{time}$, and the system becomes the following form:

$$\begin{array}{cc}\hfill D_{\ast }^{\mathsf{\alpha }}{S}^{\left(\mathsf{\alpha }\right)}& ={\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}-{\mathsf{\beta }}_{\ast }^{\mathsf{\alpha }}{S}^{\left(\mathsf{\alpha }\right)}{I}^{\left(\mathsf{\alpha }\right)}-{\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}{S}^{\left(\mathsf{\alpha }\right)}\hfill \\ \hfill D_{\ast }^{\mathsf{\alpha }}{E}^{\left(\mathsf{\alpha }\right)}& ={\mathsf{\beta }}_{\ast }^{\mathsf{\alpha }}{S}^{\left(\mathsf{\alpha }\right)}{I}^{\left(\mathsf{\alpha }\right)}-({\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}+{\mathsf{\delta }}_{\ast }^{\mathsf{\alpha }})E^{\left(\mathsf{\alpha }\right)}\hfill \\ \hfill D_{\ast }^{\mathsf{\alpha }}{I}^{\left(\mathsf{\alpha }\right)}& ={\mathsf{\delta }}_{\ast }^{\mathsf{\alpha }}{E}^{\left(\mathsf{\alpha }\right)}-({\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}+{\mathsf{\sigma }}_{\ast }^{\mathsf{\alpha }})I^{\left(\mathsf{\alpha }\right)}\hfill \\ \hfill D_{\ast }^{\mathsf{\alpha }}{R}^{\left(\mathsf{\alpha }\right)}& ={\mathsf{\sigma }}_{\ast }^{\mathsf{\alpha }}{I}^{\left(\mathsf{\alpha }\right)}-{\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}{R}^{\left(\mathsf{\alpha }\right)}.\hfill \end{array}$$

(14)

By replacing $\theta$ by ${\theta }_{\ast }^{\mathsf{\alpha }}$, the system of integer order differential equations with its epidemiological parameters follows directly from the system of fractional order differential equations when $\mathsf{\alpha }=1$.

2.4. 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 $\mathsf{\alpha }=1$ in Equation (14).

$$\begin{array}{cc}\hfill DS& =\mathsf{\mu }-\mathsf{\beta }SI-\mathsf{\mu }S\hfill \\ \hfill DE& =\mathsf{\beta }SI-(\mathsf{\mu }+\mathsf{\delta })E\hfill \\ \hfill DI& =\mathsf{\delta }E-(\mathsf{\mu }+\mathsf{\sigma })I\hfill \\ \hfill DR& =\mathsf{\sigma }I-\mathsf{\mu }R\hfill \end{array},$$

(15)

where $\mathsf{\mu },\mathsf{\beta },\mathsf{\delta },$ and $\mathsf{\sigma }$ are the model parameters described above. They all have dimensions given by $\frac{1}{\mathrm{time}}$. The last Equation in (15) is redundant since $R=1-(S+E+I)$.

2.5. Measles Model via $\mathsf{\alpha }$-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 $\mathsf{\alpha }=1$; however, in the FDE case, $\mathsf{\alpha }$ 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 $\mathsf{\alpha }$ in the parameter values. We call this new system the $\mathsf{\alpha }$-dependent ODE. By dropping the $\mathsf{\alpha }$ order derivative from the left side and $\mathsf{\alpha }$ power from $S^{\left(\mathsf{\alpha }\right)}$, $E^{\left(\mathsf{\alpha }\right)}$, and $I^{\left(\mathsf{\alpha }\right)}$ of Equation (14), our $\mathsf{\alpha }$-dependent ODE takes the following form:

$$\begin{array}{cc}\hfill DS& ={\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}-{\mathsf{\beta }}_{\ast }^{\mathsf{\alpha }}SI-{\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}S\hfill \\ \hfill DE& ={\mathsf{\beta }}_{\ast }^{\mathsf{\alpha }}SI-({\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}+{\mathsf{\delta }}_{\ast }^{\mathsf{\alpha }})E\hfill \\ \hfill DI& ={\mathsf{\delta }}_{\ast }^{\mathsf{\alpha }}E-({\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}+{\mathsf{\sigma }}_{\ast }^{\mathsf{\alpha }})I\hfill \\ \hfill DR& ={\mathsf{\sigma }}_{\ast }^{\mathsf{\alpha }}I-{\mathsf{\mu }}_{\ast }^{\mathsf{\alpha }}R.\hfill \end{array}$$

(16)

Here, the above systems may differ from the classical ODE (Equation (15)) if $\mathsf{\alpha }\ne 1$. In this case, the $\mathsf{\alpha }$ value used in the above $\mathsf{\alpha }$-dependent ODE are obtained from data fitting procedures using the FDE model (Equation (14)).

2.6. 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 $\mathsf{\alpha }$ from $S^{\left(\mathsf{\alpha }\right)}$, $E^{\left(\mathsf{\alpha }\right)}$, and $I^{\left(\mathsf{\alpha }\right)}$, for brevity.

The following two lemmas of asymptotic behavior of FDEs are given here and their proof in Appendix A for completeness.

Lemma 1.

The closed simplex region $M=\{(S,E,I)\in {\mathbb{R}}_{+}^3:0\le S+E+I\le 1\}$ is a positive invariant set for the FDE model in Equation (14).

We can find the model’s equilibrium points by setting $D_{\ast }^{\mathsf{\alpha }}S=0$, $D_{\ast }^{\mathsf{\alpha }}E=0$, and $D_{\ast }^{\mathsf{\alpha }}I=0$. Thus, there are two equilibria to the measles SEIR model (14). They are:

• the disease free equilibrium $DFE\equiv (1,0,0)$ and
• the endemic equilibrium

$$EE=(s^{\ast },e^{\ast },i^{\ast })\equiv \left(\frac{1}{R_0},\frac{\mathsf{\mu }}{\mathsf{\delta }+\mathsf{\mu }}\left(1-\frac{1}{R_0}\right),\frac{\mathsf{\mu }}{\mathsf{\beta }}(R_0-1)\right),$$

(17)

where the basic reproduction number is $R_0:=\frac{\mathsf{\beta }\mathsf{\delta }}{(\mathsf{\mu }+\mathsf{\sigma })(\mathsf{\mu }+\mathsf{\delta })}$. $EE$ exists only when $1<R_0<1+\frac{\mathsf{\beta }}{\mathsf{\mu }}$.

Remark 1.

An equilibrium is locally asymptotically stable if the eigenvalues of the Jacobian matrix of the n-dimensional system, namely ${\mathsf{\lambda }}_1,{\mathsf{\lambda }}_2,\dots ,{\mathsf{\lambda }}_n$, have the property that $|arg\left({\mathsf{\lambda }}_i\right)|>\frac{\mathsf{\alpha }\pi }{2}$, for $i=1,2,\dots ,n$, [25] (p. 158).

Thus, in general, the stability of the ordinary differential equations model implies stability of its fractional counter model of order $\mathsf{\alpha }$. 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 ${\mathsf{\lambda }}_{\ell }$ and $u_{\ell }$ for $\ell =1,2,\dots ,n$ be the eigenvalues and their respective eigenvectors of an $n\times n$ matrix A. The general solution of initial value problem consisting of a system of n linear fractional differential equations $D_{\ast }^{\mathsf{\alpha }}x\left(t\right)=Ax\left(t\right)$, such that $x\left(0\right)=x_0$, can be found to be

$$x\left(t\right)=\sum _{\ell =1}^n{c}_{\ell }{u}_{\ell }{E}_{\mathsf{\alpha }}\left({\mathsf{\lambda }}_{\ell }{t}^{\mathsf{\alpha }}\right)$$

(18)

for certain constants $c_{\ell }\in \mathbb{C}$ for $\ell =1,2,\dots ,n$, such that ${\sum }_{\ell =1}^n{c}_{\ell }{u}_{\ell }=x_0$,25] (Theorem 7.13). In case that $\mathsf{\alpha }=1$, we recover the known solution of the system of ODEs given by

$$x\left(t\right)=\sum _{\ell =1}^n{c}_{\ell }{u}_{\ell }exp\left({\mathsf{\lambda }}_{\ell }t\right).$$

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 ${\mathsf{\lambda }}_2$ and ${\mathsf{\lambda }}_3$, are conjugate complex-valued. In that situation, $x\left(t\right)$ would oscillate with inter-peak periods, called inter-epidemic period in disease modeling, given by $2\pi {\left(Im\left({\mathsf{\lambda }}_2\right)\right)}^{-1}$ [14]. If $Re\left({\mathsf{\lambda }}_{\ell }\right)<0$ for all ℓ, then the oscillations will be damped to zero. That damped oscillation is straightforward in the case of $\mathsf{\alpha }=1$ due to the exponential damping in the superposition of the sine and cosine functions. That behavior, however, is not straightforward for $0<\mathsf{\alpha }<1$.

2.7. Numerical Simulations

Since the mean of the subordinator process is $\mathbb{E}\left({\mathcal{T}}_{\mathsf{\alpha }}\left(t\right)\right)={\displaystyle \frac{t^{\mathsf{\alpha }}}{\mathrm{\Gamma }(\mathsf{\alpha }+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:

$$\begin{array}{cc}\hfill D_{\ast }^{\mathsf{\alpha }}x\left(t\right)& =f(t,x(t\left)\right),\mathrm{for}0<t\le T,\hfill \\ \hfill x\left(0\right)& =x_0,\hfill \end{array}$$

(19)

for some $T>0$. A solution of Equation (19) is approximated by the deterministic time subordination

$$x\left(t\right)=y\left(\frac{t^{\mathsf{\alpha }}}{\mathrm{\Gamma }(\mathsf{\alpha }+1)}\right),$$

(20)

of $y\left(s\right)$, the solution of the ordinary differential equation

$$\begin{array}{cc}\hfill \frac{dy\left(s\right)}{ds}& =g(s,y\left(s\right)),\mathrm{for}0<s\le \frac{t^{\mathsf{\alpha }}}{\mathrm{\Gamma }(\mathsf{\alpha }+1)}\hfill \\ \hfill y\left(0\right)& =x_0.\hfill \end{array},$$

(21)

where

$$\begin{array}{c}\hfill g(s,y\left(s\right))=f(t-{(t^{\mathsf{\alpha }}-s\mathrm{\Gamma }(\mathsf{\alpha }+1))}^{\frac{1}{\mathsf{\alpha }}},x(t-{(t^{\mathsf{\alpha }}-s\mathrm{\Gamma }(\mathsf{\alpha }+1))}^{\frac{1}{\mathsf{\alpha }}}))\end{array}$$

(22)

for all $0<t\le 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.

Algorithm 1 Numerical solution of $D_{\ast }^{\mathsf{\alpha }}x\left(t\right)=f(t,x\left(t\right))$ for $0<t<T$ with $x\left(0\right)=x_0$.

Input: $\mathsf{\alpha },T,f(t,x(t\left)\right),m,n$ Output: $x\left(t\right)$

begin

Divide the interval $[0,T]$ into n sub-intervals using $$0=t_0<t_1<\dots <t_n=T.$$

for$i=1,2,\dots ,n$

Divide the interval $[0,{\displaystyle \frac{t_i^{\mathsf{\alpha }}}{\mathrm{\Gamma }(\mathsf{\alpha }+1)}}]$ into further m sub-intervals using $$0=s_0<s_1<\dots <s_m={\displaystyle \frac{t_i^{\mathsf{\alpha }}}{\mathrm{\Gamma }(\mathsf{\alpha }+1)}}.$$

Solve the system $Dy\left(s\right)=f(t_i-{(t_i^{\mathsf{\alpha }}-s\mathrm{\Gamma }(\mathsf{\alpha }+1))}^{\frac{1}{\mathsf{\alpha }}},y\left(s\right))$ with $y\left(0\right)=x_0$ using Euler or Runge-Kutta methods on $s_0,s_1,\dots ,s_m$.

Retain $x\left(t_i\right)=y\left(s_m\right)$.

end

Return $[x_0,x\left(t_1\right),x\left(t_2\right),\dots ,x\left(t_n\right)]$.

end

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

$$L(\mathsf{\alpha },\mathsf{\beta },\mathsf{\mu },\mathsf{\delta },\mathsf{\sigma })=\sum _{i=1}^n{(I_i-{\widehat{I}}_i)}^2$$

for $\mathsf{\alpha }\in (0,1]$, and $\mathsf{\beta },\mathsf{\mu },\mathsf{\delta },\mathsf{\sigma }\in (0,\infty )$, where $\{I_i,i=1,\dots ,n\}$ is the data of actual proportion of infections and $\{{\widehat{I}}_i,i=1,\dots ,n\}$ is the simulated proportion of infections. The values ${\widehat{I}}_i$ approximating $I\left(t_i\right)$ 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.

3. 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 $\mathsf{\alpha }$. For smaller $\mathsf{\alpha }$, 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 $\mathsf{\alpha }\to 1$. To further compare the two modeling approaches, we consider the analogue ODEs derived for specific $\mathsf{\alpha }$ 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.

The estimates of $\mathsf{\alpha }$ are $0.999999999524,0.999999999614,$ and $0.875784505546$ for New York, London, and Portsmouth, respectively. The values obtained for $\mathsf{\alpha }$ for the New York and London data sets are very close to one, where the FDE model converges to the ODE model (Equation (15)). Statistical measurements on the data fittings are presented in

Table A4. Parameter estimation and statistical measures of fittings for the classical ODE model and FDE model using different data sets.

Data	Model	Estimated Parameters, $(\mathsf{\alpha },\mathsf{\mu },\mathsf{\beta },\mathsf{\delta },\mathsf{\sigma })$	SSE	AIC	BIC
New York	ODE	$(NA,0.0028,119.22,16.73,10.19)$	$1.29\times {10}^6$	250.54	253.87
	FDE	$(0.999999999524,0.0029,116.34,19.39,10.37)$	$1.34\times {10}^6$	255.36	259.53
London	ODE	$(NA,6.79\times {10}^{-4},153.44,1.99,4.27)$	$2.01\times {10}^6$	389.36	394.54
	FDE	$(0.999999999614,0.0015,122.6849,2.18,8.75)$	$2.04\times {10}^6$	393.34	398.85
Portsmouth	ODE	$(NA,{10}^{-6},228.61,0.46,3.33)$	$4.57\times {10}^4$	277.94	282.98
	FDE	$(0.875784505546,2.52\times {10}^{-4},278.72,1.52,5.24)$	$3.22\times {10}^4$	271.92	278.21

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

4. 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 $\mathsf{\alpha }$ 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 $\mathsf{\alpha }=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 ${\mathcal{T}}_{2\mathsf{\alpha }}$, 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.

5. 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, $\mathsf{\alpha }=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.