Global Stability of an SEIR Epidemic Model where Empirical Distribution of Incubation Period has Approximated by Coxian Distribution

TL;DR

This paper introduces a Coxian distributed SEIR epidemic model that incorporates empirical incubation period data, demonstrating that its global dynamics depend solely on the basic reproduction number and offering a flexible approach for realistic epidemic modeling.

Contribution

The paper develops a Coxian distributed SEIR model that accurately incorporates empirical incubation periods, enhancing classical models' realism and analytical tractability.

Findings

Global dynamics depend on the basic reproduction number.

Application to empirical data demonstrates model's practical utility.

Coxian distribution's flexibility allows approximation of various incubation period distributions.

Abstract

In this work, we have developed a Coxian distributed SEIR model in incorporating an empirical incubation period. We show that the global dynamics are completely determined by a basic reproduction number. An application of the Coxian distributed SEIR model using data of an empirical incubation period is explored. The model may be useful for resolving causing the realistic intrinsic parts in classical epidemic models since Coxian distribution approximately converges to any distribution.

Tables

Table 1. Table 1: SEIR models generated by various distributions with delay effect in the previous papers.

Distribution, Survival function($P_{△}\left(t\right)$)	Derived model ($S\left(t\right)$, $I\left(t\right)$)	Basic reproduction number ($\beta \left(t\right)\equiv \overline{\beta }$)	Characteristics
Dirac-delta $P_D\left(t\right)=\{\begin{array}{cc}1,u\in [0,\tau ),& \\ 0,u\in [\tau ,\infty ),& \end{array}$ for $\tau >0$, constant	$\frac{\text{d}E\left(t\right)}{\text{d}t}=\beta \left(t\right)S\left(t\right)I\left(t\right)-\beta \left(t-\tau \right)S\left(t-\tau \right)I\left(t-\tau \right)\exp \left(-\mu \tau \right)-\mu E\left(t\right),$ $\frac{\text{d}I\left(t\right)}{\text{d}t}=\beta \left(t-\tau \right)S\left(t-\tau \right)I\left(t-\tau \right)\exp \left(-\mu \tau \right)-\left(\gamma +\mu \right)I\left(t\right).$	$\frac{\overline{\beta }\exp \left(-\mu \tau \right)}{\gamma +\mu }$	System of DDE, Exact time of "delay", NOT distributed
Exponential $P_E\left(t\right)=\exp \left(-\lambda t\right)$	$\frac{\text{d}E\left(t\right)}{\text{d}t}=\beta \left(t\right)S\left(t\right)I\left(t\right)-\left(\lambda +\mu \right)E\left(t\right),$ $\frac{\text{d}I\left(t\right)}{\text{d}t}=\lambda E\left(t\right)-\left(\gamma +\mu \right)I\left(t\right).$	$\frac{\overline{\beta }\mu }{\left(\lambda +\mu \right)\left(\gamma +\mu \right)}$	System of simple ODE, Memoryless property
Gamma $P_G\left(u\right)={\sum }_{i=0}^{n-1}\frac{1}{i!}{e}^{-n\lambda u}{\left(n\lambda u\right)}^i$, for $n$, positive integer	$\frac{\text{d}{E}_n\left(t\right)}{\text{d}t}=\beta \left(t\right)S\left(t\right)I\left(t\right)-\left(n\lambda +\mu \right)E_n,$ $\frac{\text{d}{E}_i\left(t\right)}{\text{d}t}=n\lambda E_{i+1}-\left(n\lambda +\mu \right)E_i,\text{ for }i=n-1,\cdots ,2,1,$ $\frac{\text{d}I\left(t\right)}{\text{d}t}=n\lambda E_1\left(t\right)-\left(\gamma +\mu \right)I\left(t\right),$ with $E\left(t\right)={\sum }_{i=0}^{n-1}{E}_{n-i}\left(t\right)$.	$\frac{\overline{\beta }}{\gamma +\mu }{\left(\frac{n\lambda }{n\lambda +\mu }\right)}^n$	Linear chain, unimodal, short-tailed distribution
Mittag-Leffler $P_M\left(t\right)=E_{\alpha ,1}\left(-{\left(t/\zeta \right)}^{\alpha }\right),$ for where $E_{\alpha ,\beta }\left(z\right)={\sum }_{k=0}^{\infty }\frac{z^k}{\Gamma \left(\alpha k+\beta \right)}$	$\frac{\text{d}E\left(t\right)}{\text{d}t}=\beta \left(t\right)S\left(t\right)I\left(t\right)-\exp \left(-\mu t\right){\zeta }^{-\alpha }{}_0{}^{}𝒟_t^{1-\alpha }\left[\exp \left(\mu t\right)E\left(t\right)\right]-\mu E\left(t\right),$ $\frac{\text{d}I\left(t\right)}{\text{d}t}=\exp \left(-\mu t\right){\zeta }^{-\alpha }{}_0{}^{}𝒟_t^{1-\alpha }\left[\exp \left(\mu t\right)E\left(t\right)\right]-\left(\gamma +\mu \right)I\left(t\right),$ where ${}_0{}^{}𝒟_t^{1-\alpha }\left[f\left(t\right)\right]=\frac{1}{\Gamma \left(\alpha \right)}\frac{\text{d}}{\text{d}t}{\int }_0^t{\left(t-u\right)}^{\alpha -1}f\left(u\right)\text{d}u$.	$\frac{\overline{\beta }}{\gamma +\mu }\cdot \frac{1}{1+{\left(\zeta \mu \right)}^{\alpha }}$	System of FDE, heavy-tailed distribution, Hard to get exact form of distribution

Equations

$$\frac{\text{d}S(t)}{\text{d}t}$$

$$E(t)$$

$$\frac{\text{d}I(t)}{\text{d}t}$$

$$\frac{\text{d}R(t)}{\text{d}t}$$

$$P(u)=\mathbb{P}(T>u)=1-F(u),$$

$$\begin{array}{rl}\mathbf{p}& =\left[\begin{array}{cccc}1& 0& \cdots & 0\end{array}\right]\in M_{1\times n}(\mathbb{R}),\\ \mathbf{Q}& =\left[\begin{array}{cccccc}-{\lambda }_n& {\bar{p}}_{n-1}{\lambda }_n& 0& & \cdots & 0\\ 0& -{\lambda }_{n-1}& {\bar{p}}_{n-2}{\lambda }_{n-1}& 0& \cdots & 0\\ & & & ⋮& & \\ 0& \cdots & & & -{\lambda }_2& {\bar{p}}_1{\lambda }_2\\ 0& \cdots & & & 0& -{\lambda }_1\end{array}\right]\in M_{n\times n}(\mathbb{R}),\\ \mathbf{q}& =-Q1={\left[\begin{array}{ccccc}{p}_{n-1}{\lambda }_n& p_{n-2}{\lambda }_{n-1}& \cdots & p_1{\lambda }_2& {\lambda }_1\end{array}\right]}^{\top }\in M_{n\times 1}(\mathbb{R}),\end{array}$$

$$P(t)=\mathbf{p}\exp (t\mathbf{Q})1.$$

$$[\mathcal{L}{f}_X](s)=\mathbf{p}(s\mathbf{I}-\mathbf{Q}{)}^{-1}\mathbf{q}=\sum _{i=0}^{n-1}\left[a_i\prod _{j=1}^{n-i}\frac{{\lambda }_{n+1-j}}{s+{\lambda }_{n+1-j}}\right],$$

$$\mathbf{P}(t)$$

$$P(t)=\sum _{i=1}^n{P}_n(t).$$

$$\begin{array}{rl}\frac{\text{d}{P}_n}{\text{d}t}& =-{\lambda }_n{P}_n,\\ \frac{\text{d}{P}_i}{\text{d}t}& ={\bar{p}}_i{\lambda }_{i+1}{P}_{i+1}-{\lambda }_i{P}_i,\text{ for }i=n-1,\cdots ,1,\end{array}$$

$$\frac{\text{d}P(t)}{\text{d}t}=\sum _{i=1}^n\frac{\text{d}{P}_i(t)}{\text{d}t}=-{\lambda }_1{P}_1(t)-\sum _{i=1}^{n-1}{p}_i{\lambda }_{i+1}{P}_{i+1}(t).$$

$$E(t)=\sum _{i=1}^n{\int }_0^{\infty }{P}_i(u)\beta (t-u)S(t-u)I(t-u)\exp (-\mu u)\text{d}u.$$

$$E_i(t)={\int }_0^{\infty }{P}_i(u)\beta (t-u)S(t-u)I(t-u)\exp (-\mu u)\text{d}u,$$

$$\begin{array}{rl}\frac{\text{d}{E}_n(t)}{\text{d}t}& =\beta (t)S(t)I(t)-({\lambda }_n+\mu )E_n(t),\\ \frac{\text{d}{E}_i(t)}{\text{d}t}& ={\bar{p}}_i{\lambda }_{i+1}{E}_{i+1}(t)-({\lambda }_i+\mu )E_i(t),i=n-1,\cdots ,2,1,\end{array}$$

$$\frac{\text{d}I(t)}{\text{d}t}={\int }_0^t\beta (t-u)S(t-u)I(t-u)\left[{\lambda }_1{P}_1(u)+\sum _{i=1}^{n-1}{p}_i{\lambda }_{i+1}{P}_{i+1}(u)\right]\exp (-\mu u)\text{d}u-(\gamma +\mu )I(t),$$

$$\frac{\text{d}I(t)}{\text{d}t}={\lambda }_1{E}_1(t)+\sum _{i=1}^{n-1}{p}_i{\lambda }_{i+1}{E}_{i+1}(t)-(\gamma +\mu )I(t).$$

$$\begin{array}{rl}\frac{\text{d}S(t)}{\text{d}t}& =\mu -\beta (t)S(t)I(t)-\mu S(t),\\ \frac{\text{d}{E}_n(t)}{\text{d}t}& =\beta (t)S(t)I(t)-({\lambda }_n+\mu )E_n(t),\\ \frac{\text{d}{E}_i(t)}{\text{d}t}& ={\bar{p}}_i{\lambda }_{i+1}{E}_{i+1}(t)-({\lambda }_i+\mu )E_i(t),i=n-1,\cdots ,2,1,\\ \frac{\text{d}I(t)}{\text{d}t}& =\sum _{i=0}^{n-1}{p}_i{\lambda }_{i+1}{E}_{i+1}(t)-(\gamma +\mu )I(t),\\ \frac{\text{d}R(t)}{\text{d}t}& =\gamma I(t)-\mu R(t),\end{array}$$

$$\Omega =\left\{(S,E_n,E_{n-1},\cdots ,E_1,I,R)\in {\mathbb{R}}_{+}^{n+3}:S+\sum _{j=1}^n{E}_j+I+R\equiv 1\right\}.$$

$$\frac{\text{d}S(t)}{\text{d}t}=\mu -(\zeta (t)+\mu )S(t)\ge -(\zeta (t)+\mu )S(t)$$

$$S(t_1)\ge S(0)\exp \left\{-\mu t_1-{\int }_0^{t_1}\zeta (u)\text{d}u\right\}>0$$

$$\hat{P}={\int }_0^{\infty }\exp (-\mu u)P(u)\text{d}u.$$

$$0\le \hat{P}\le {\int }_0^{\infty }\exp (-\mu u)\text{d}u=1/\mu$$

$$R_0=\bar{\beta }\cdot (1-\mu \hat{P})\cdot \frac{1}{\gamma +\mu }.$$

$$1-\mu \hat{P}={\int }_0^{\infty }\exp (-\mu u)f_X(u)\text{d}u=[\mathcal{L}{f}_X](\mu ).$$

$$R_0=\frac{\bar{\beta }}{\gamma +\mu }\sum _{i=0}^{n-1}\left[a_i\prod _{k=1}^{n-i}\frac{{\lambda }_{n+1-k}}{{\lambda }_{n+1-k}+\mu }\right],$$

$${\mathcal{E}}^d=(S^d,E_n^d,E_{n-1}^d,\cdots ,E_1^d,I^d,R^d)=(1,0,0,\cdots ,0),$$

$${\Omega }_d=\left\{(S,E_n,E_{n-1},\cdots ,E_1,I,R)\in {\mathbb{R}}_{+}^{n+3}:E_n=E_{n-1}=\cdots =E_1=I=R=0\right\}.$$

$$V=\sum _{j=1}^n{w}_j{E}_j+I,$$

$$w_j=\sum _{l=1}^j{\alpha }_l^j\prod _{k=l}^j\frac{{\lambda }_k}{{\lambda }_k+\mu },$$

$${\alpha }_l^j=\{\begin{array}{ll}{p}_{j-1}& \text{if l=j}\\ p_{l-1}{\prod }_{\xi =l}^{j-1}{\bar{p}}_{\xi }& \text{otherwise, i.e. }l=1,2,\cdots ,j-1.\end{array}$$

$$\begin{array}{rl}\frac{\text{d}V}{\text{d}t}& =\sum _{j=1}^n{w}_j\frac{\text{d}{E}_j}{\text{d}t}+\frac{\text{d}I}{\text{d}t},\\ & =w_n\left[\bar{\beta }S(t)I(t)-({\lambda }_n+\mu )E_n(t)\right]+\sum _{j=1}^{n-1}{w}_j\left[{\bar{p}}_j{\lambda }_{j+1}{E}_{j+1}(t)-({\lambda }_j+\mu )E_j(t)\right]\\ & +\sum _{j=0}^{n-1}{p}_j{\lambda }_{j+1}{E}_{j+1}(t)-(\gamma +\mu )I(t),\\ & =w_n\bar{\beta }S(t)I(t)-(\gamma +\mu )I(t)+\sum _{j=2}^n\left[-w_j({\lambda }_j+\mu )+w_{j-1}{\bar{p}}_{j-1}{\lambda }_j+p_{j-1}{\lambda }_j\right]E_j(t)\\ & +[p_0{\lambda }_1-w_1({\lambda }_1+\mu )]E_1(t),\end{array}$$

Full text

K****eywords Basic reproduction number  $\cdot$ Coxian distributed SEIR model  $\cdot$ Global stability $\cdot$ Infectious diseases modeling  $\cdot$ LaSalle’s invariance principle

1 Introduction

Compartmental models in epidemiology are well used to simplify mathematical modeling of infectious diseases ([1]). Its origin is in the early 20th century, with an important early work being that of Kermack and McKendrick in 1927 ([2]). The population is subdivided into three groups: susceptible, infectious and recovered groups, and the model is referred to as an SIR epidemic model. The model is applicable to infectious diseases such as measles, chickenpox, or mumps, where infection confers immunity.

For many important infections there is a significant incubation period during which the individual has been infected but is not yet infectious themselves ([3]). During this period the individual is in compartment $E$, known as the exposed term, and such type models are called SEIR models. To express the distribution of incubation period, we consider age-structured SEIR framework. Let $S(t)$, $E(t)$, $I(t)$ and $R(t)$ be the fraction of susceptible, exposed, infectious and recovered population at time $t$, respectively. Motivated by [4], consider an SEIR model of integro-differential type:

[TABLE]

where $\beta(t)$ is the time-dependent contact rate, $1/\gamma$ is the mean period of infectiousness, $P(u)$ is the probability of remaining incubated (assumed asymptomatic and not infectious) to time $u$ from entering, and $\mu$ is the natural mortality. $P_{u}$ denotes the derivative with respect to $u$. The initial conditions are given as $S(0)=S_{0},~{}E(0)=E_{0},~{}I(0)=I_{0},~{}R(0)=R_{0}\textrm{, and }S_{0}+E_{0}+I_{0}+R_{0}=1$. In general, it is reasonable to assume that a function $P:[0,\infty)\rightarrow[0,1]\in\mathcal{L}^{1}$ is Lebesgue integrable, differentiable a.e., non-increasing, $P(0)=1$, and $\lim_{t\rightarrow\infty}P(t)=0$ ([5]). Especially, let $T$ be a continuous random variable of incubation period of some infectious diseases. If we know the distribution of incubation period heuristically, we can define $P$ by

[TABLE]

where $F$ is a cumulative density function of random variable of infectious period $T$. $P$ is called the survival function. Many previous studies have developed with the various survival functions $P$ such as Dirac-delta [6, 7], Exponential [8], Gamma [9, 10, 11, 12], Mittag-Leffler [13, 14], and joint [15, 16] types. In some cases, this framework of SEIR model can have been simplified by system of ordinary or fractional differential equations with various incubation period distribution (see Table 1).

Following these, in this work, we construct an SEIR epidemic model, which is named Coxian distributed SEIR model, to incorporate empirical incubation period distribution. Coxian distribution, one of the phase-type distributions, can be considered as a mixture of hypo– and hyper–exponential distributions. The novelty of the distribution is the density in the class of all non–negative distribution functions [17] and so all types of incubation periods are approximated to Coxian distribution. A basic reproduction number, $R_{0}$, was derived. Applying a Lyapunov theory, we show that the basic reproduction number determines the global stability of equilibrium of our model with constant contact rate case. We also give an application illustrating how to use a Coxian distributed SEIR model with empirical incubation period data. The limitation of the model is also discussed in the last section.

2 Derivation of Coxian Distributed SEIR Model

A Coxian distribution with $n$ phases is defined as the time until absorption into state [math], starting from state $n$, of the Markov process with discrete states in continuous time. A focal state, $X$, can be partitioned into $n$ sub-states $X_{n},\cdots,X_{1}$, each with independent dwell time distributions that are exponentially distributed with rates $\lambda_{i}$, $i=n,n-1,\cdots,1$. Inflow rates into the state $X_{n}$ can be described by non-negative, integrable inflow rates ${I}_{X}(t)$. Particles that transition out of a sub-state $X_{k}$ at time $t$ transition into either a different sub-state $X_{k-1}$ with probability $\bar{p}_{k-1}$, or enter the recipient state $X_{0}$ with probability $p_{k-1}=1-\bar{p}_{k-1}$, for $k=n,n-1,\cdots,1$. Then, Coxian distributed random variable $X$ [17] is the time until absorption in $X_{0}$ starting from state $n$ and the probability density function $f_{X}$ is given as $f_{X}(t)=\mathbf{p}\exp{(t\mathbf{Q})}\mathbf{q}$ where

[TABLE]

$\mathbf{1}$ is an $n\times 1$ vector of ones and $\mathbf{Q}$ is transition rate matrix, $0<p_{i}<1$ for $i=1,2,\cdots,n-1$ and $\ \lambda_{i}>0$ for $i=1,2,\cdots,n$. The survival function $P$ is given as

[TABLE]

The Laplace transform $\mathcal{L}$ of $f_{X}$ is given explicitly by

[TABLE]

where $\mathbf{I}$ is $n\times n$ identity matrix and $a_{i}=p_{i}\prod_{j=1}^{n-1-i}\bar{p}_{j+i}$.

Derivation of Coxian distributed SEIR model is the following: assume that infectious period is distributed by Coxian with the survival function (4) in the model (12). First, note that if we put the $1\times n$ vector $\mathbf{P}$ as:

[TABLE]

where $P_{i}$’s are functions of $t$, then the survival function $P$ is expressed as

[TABLE]

Since $\frac{\textrm{d}[\exp{(t\mathbf{Q})}]}{\textrm{d}t}=\exp{(t\mathbf{Q})}\mathbf{Q}$, $\mathbf{P}^{\prime}(t)=\mathbf{P}(t)\mathbf{Q}$ holds and so

[TABLE]

and

[TABLE]

Now, if we put (4) to (1b), (1b) has the form from (6):

[TABLE]

So if we put

[TABLE]

for $i=n,n-1,\cdots,1$ and differentiating (9) for each $i=n,n-1,\cdots,1$, then (7) yields:

[TABLE]

and substituting (8) into (1c) yields

[TABLE]

and so the (1c) becomes

[TABLE]

From (10) and (11), targeted $n-$chained Coxian distributed SEIR model is derived from (1) when we put $p_{0}=1$, as follows:

[TABLE]

with $\sum_{i=1}^{n}E_{i}(t)=E(t)$. The schematic diagram for the model (12) is depicted in Figure 1.

Basic properties of the model are below: the region

[TABLE]

is obviously positive-invariant. Also if the initial data $S(0),E_{i}(0),I(0)$ and $R(0)$ for $i=n,n-1,\cdots,1$ are positive, then the solutions $S(t),E_{i}(t),I(t)$ and $R(t)$ of the model (12) are non-negative for all $t>0$. Indeed, if we assume that there exists $t_{1}>0$ such that $S(t)>0,E_{i}(t)>0,I(t)>0,R(t)>0$, $i=n,n-1,\cdots,1$ in $t\in[0,t_{1})$ and $S(t_{1})\cdot\prod_{i=1}^{n}E_{i}(t_{1})\cdot I(t_{1})\cdot R(t_{1})=0$ (i.e., one of states is 0 at the time $t_{1}$). Then,

[TABLE]

when we take a force of infection $\zeta(t)=\beta(t)I(t)$, we may get

[TABLE]

holds from Gronwall. Similarly, it can be shown that $E_{i}(t_{1})>0,I(t_{1})>0$ and $R(t_{1})>0$ for $i=n,n-1,\cdots,1$ and this yields a contradiction.

3 The Basic Reproduction Number, $R_{0}$

When the parameters are constant, we can compute the basic reproduction number. Assume the number of infectious human is rare in the early phase, i.e. $S\approx 1$. Consider $\beta(t)\equiv\bar{\beta}$. Define the average time that an individual remains in the exposed class becoming infectious without dying, $\hat{P}$, as

[TABLE]

We note that $\mu\hat{P}$ represents the probability that an exposed individual will die during the course of incubating, and so the probability of surviving the exposed class is $1-\mu\hat{P}$. Note that since $0\leq P(t)\leq 1$ for $t\geq 0$, the inequality

[TABLE]

holds and so the probability $\mu\hat{P}$ is well-defined. From the definition of the basic reproduction number, we can represent the basic reproduction number, denoted by $R_{0}$, as

[TABLE]

Since $\mu\hat{P}=1+\int_{0}^{\infty}\exp{(-\mu u)}P^{\prime}(u)\textrm{d}u$ from integration by parts, and the relation $P^{\prime}(t)=-f_{X}(t)$ in the usual sense, we get

[TABLE]

From (5), $R_{0}$ is explicitly represented as:

[TABLE]

where $a_{i}=p_{i}\prod_{j=1}^{n-1-i}\bar{p}_{j+i}$. Note that if we ignore demographic part, i.e. $\mu=0$, then $R_{0}$ is expressed as $\bar{\beta}/\gamma$, same as the one in classical model.

4 Equilibria Analysis for the case of constant contact rate

Consider the case when $\beta(t)=\bar{\beta}$, a positive constant. The model (12) has a equilibrium

[TABLE]

which is called the disease free equilibrium. For convenience, let

[TABLE]

Now we prove the following properties:

Theorem 1**.**

The disease free equilibrium, $\mathcal{E}^{d}$, of the model (12) is globally asymptotically stable in the domain $\Omega_{d}$ if $R_{0}\leq 1$.

Proof.

Consider the function

[TABLE]

for positive constants $w_{j}$, for $j=1,\cdots,n$, where

[TABLE]

with

[TABLE]

for fixed $j$. Then, we can note that $w_{n}={(\gamma+\mu)R_{0}}/{\bar{\beta}}$. The derivative of $V$ to time $t$ is given by

[TABLE]

by rearrangement. Moreover, since

[TABLE]

and

[TABLE]

and so we can get the recurrence relation:

[TABLE]

for each $j=2,3,\cdots,n$. Thus, if we substitute (16) and (17) into (15), we may get

[TABLE]

Therefore, $V^{\prime}\leq 0$ if $R_{0}\leq 1$ and $V^{\prime}=0$ if and only if $E_{n}=E_{n-1}=\cdots=E_{1}=I=0$. Also, substituting $I=0$ into the equations for $S$ and $R$ shows that $S\rightarrow 1$ and $R\rightarrow 0$ as $t\rightarrow\infty$. Hence, $V$ is the Lyapunov function in the domain $\Omega_{d}$. From the LaSalle’s invariance principle, the disease free equilibrium $\mathcal{E}^{d}$ is globally asymptotically stable in $\Omega_{d}$ if $R_{0}\leq 1$. ∎

Next we consider the endemic equilibrium, $\mathcal{E}^{\ast}=(S^{\ast},E_{n}^{\ast},\cdots,E_{1}^{\ast},I^{\ast},R^{\ast})$, whose components are all positive. Solving the algebraic system in (12) at fixed point gives

[TABLE]

From (18a) and (18b),

[TABLE]

Substituting (19) into (18d) yields

[TABLE]

Moreover, if we substitute (20) into (18a), then

[TABLE]

since $S^{\ast}\neq 1$. Hence $I^{\ast}>0$ if $R_{0}>1$. Also, if $I^{\ast}>0$, then $0<S^{\ast}=1/R_{0}<1$ and $R^{\ast}>0$ from (18e). Moreover, $E_{i}$’s are positive since $1-S^{\ast}>0$. Conversely, if $R_{0}\leq 1$, then the model has no positive equilibrium. These results are summarized below.

Theorem 2**.**

The endemic equilibrium $\mathcal{E}^{\ast}$ of the model (12) exists uniquely when $R_{0}>1$, and no endemic equilibrium otherwise.

Finally, we claim the following:

Theorem 3**.**

The endemic equilibrium, $\mathcal{E}^{\ast}$, of the model (12) is globally asymptotically stable in the interior of $\Omega$ if $R_{0}>1$.

Proof.

For convenience, define

[TABLE]

where $f_{X}(x)=-P^{\prime}(x)$. Then, $1-\mu\hat{P}=\hat{f}(0)$. Motivated by [18], consider a Lyapunov function $V_{e}\equiv V_{e}(S,I)$ as

[TABLE]

where

[TABLE]

and $G(x)=x-1-\ln x$. Notice that $V_{e}=0$ when $S=S^{\ast},E_{n}=E_{n}^{\ast},\cdots,R=R^{\ast}$ and $V_{e}>0$ otherwise. Differentiating $V_{1e}$ and $V_{2e}$ with respect to the time $t$ yields

[TABLE]

from the relation (18a) and (18d), and,

[TABLE]

Combining (21) and (22), we can get

[TABLE]

since $G(x)\geq 0$ for all $x>0$ and so the integrand of the last term are positive. Note $V_{e}^{\prime}=0$ only if at the endemic equilibrium and $V_{e}^{\prime}<0$ otherwise. Therefore, by the LaSalle’s principle, every solution to the equation in the model (12) approaches the endemic equilibrium for $R_{0}>1$. Thus the endemic equilibrium $\mathcal{E}^{\ast}$ is globally asymptotically stable if $R_{0}>1$. ∎

Theorem 1 and Theorem 3 indicate that the disease could be eliminated by maintaining the basic reproduction number less than unity, and conversely, when $R_{0}$ is greater than unity, the disease persists in the epidemiological point of view.

5 An application

The 2009 epidemic of influenza H1N1 in Canada is investigated to explain the procedure from fitting incubation period to applying to our model. The procedure is summarized: approximating empirical distribution to Coxian distribution, and investigating a Coxian distributed SEIR model.

First, we approximate the empirical distribution of incubation period to Coxian distribution, $1-P$, as denoted in (2). Empirical data of incubation period are captured from [19] and the data are fit to the Exponential and Coxian, respectively. Unfortunately, there is no criterion for choosing the number of substates $n$. However, Akaike Information Criterion, AIC, gives the relative quality of statistical models for a given set of data [20] and so we choose $n$ that makes the AIC smallest. Next, the human case data are fit to Coxian distributed SEIR model. The report [21] is referred to obtaining daily human case data. We normalized the model (12) as multiplying the human population. Since the duration of epidemic is short related to human lifespan, we ignore the demographic effect, and so $\mu=0$ is assumed. Since the rate of infection-related death is too small to be considered, we don’t consider the death from the epidemic. Motivated by [19], the mean duration of infectiousness $1/\gamma$ is assumed as $7.1$. We consider the initial time ($t=0$) as April 14 when the epidemic started. Since there are two turning points in April 29 and June 4 in the duration of epidemic, we assume the transmission rate $\beta(t)$ consists of step function:

[TABLE]

and estimate $\beta(t)$ by fitting between the case data and the cumulative prevalence, $\int_{0}^{t}\left\{\sum_{i=0}^{n-1}p_{i}\lambda_{i+1}E_{i+1}(u)\right\}\textrm{d}u$, from the model as minimizing the sum of squared residual with the corresponding data.

Figure 2a(a) shows the result of empirical distribution of incubation period to Coxian distributions of 12-chains, which is determined by AIC(corresponding AIC$=-36.3$). We see that the fit curve explain the distribution of empirical incubation period well. Figure 2b(b) indicates the result of investigating Coxian distributed SEIR model. Well fit values of set of parameters are $\beta_{0}=0.642$, $\beta_{1}=0.252$, and $\beta_{2}=0.131$ and so the basic reproduction number (13), $R_{0}(=\beta_{0}/\gamma)$ is $4.561$. To compare the result with the classical model, we fit the empirical distribution of incubation period to exponential and investigate the classical exponential SEIR model. In this case, AIC value of distribution fit is $-33.3$. Exponential model gives $\beta_{0}=0.584$, $\beta_{1}=0.234$ and $\beta_{2}=0.133$ and so the $R_{0}$ is $4.146$. From the example, we can see that investigating Coxian distributed SEIR model might give such a nice fit results, and moreover, $R_{0}$ can vary greatly as considering the realistic distribution. This result supports the fact that common assumption for exponentially distributed incubation period always underestimate the basic reproductive number of infection from onset data and considering a realistic incubation period distribution is important to modeling epidemics ([22]).

6 Discussion and Conclusion

Many previous works have strongly emphasized that modelers should be cautious for considering the intrinsic facts to classical frameworks when the epidemic models for public health is proposed. In the study, we have derived an SEIR model based on Coxian distribution which approximates the distribution of the incubation period. In mathematical analysis, we proved that global stability of equilibria with respect to the value of $R_{0}$. The novelty of the model is that all types of distributions of the incubation period are fit with Coxian distribution and in addition, since it has realistic distribution of incubation period, it enables to describe the SEIR model of a particular type of distribution, like bimodal, that is not expressed by the conventional way, such as P. vivax malaria in temperate regions. Moreover, we may extend the Coxian distributed model not just incubation period, but infectious period. The model has some limitations even if several benefits. First, a number of the parameters are needed to fit empirical data in which they brings computation works. Secondly, loss of biological meaning could be yielded because of going out to absorption state without going through the whole chain of incubation. However, we believe that the model may be applicable for resolving the problem caused by simplicity of the models in some cases when sufficient empirical information of the incubation period is given, as shown in the previous section.

Availability of data and material

The datasets analyzed during the current study are available in [19].

Competing interests

The authors declare that they have no competing interests.

Funding

Not applicable.

Author’s contributions

IHJ conceptualized the roles and designed the studies. SK did mathematical analysis of the model and illustrated an application. SK and JHB discussed and analyzed all results. All members verified the results and wrote and reviewed the manuscript.

Figures

Tables