We consider an averaging principle for the endemic SIR model in a semi-Markov random media. Under stationary conditions of a semi-Markov media we show that the perturbed endemic SIR model converges to the classic endemic SIR model with averaged coefficients. Numerical toy examples and their interpretations are also presented for two-state Markov and semi-Markov chains. We also discuss two numerical examples involving real data: 1) Dengue Fever Disease (Indonesia and Malaysia (2009)) and 2) Cholera Outbreak in Zimbabwe (2008-2009). Novelty of the paper consists in studying of an endemic SIR model in semi-Markov random media and in implementations and interpretations of the results through numerical toy examples and discussion of numerical examples with real data.

Keywords: endemic SIR model, semi-Markov random media, averaging principle, averaged endemic SIR model, two-state Markov chain, two-state semi-Markov chain, weibull disctribution, failure rate for the disease

In the last years, deterministic and stochastic epidemic models, in both discrete and continuous time, have been studied Nasel,^1,2 Ball & Lyne³ Both model types are needed, and both have their advantages and weaknesses Hethcote,⁴ Allen & Burguin.⁵ The deterministic models lead to powerful qualitative results with important threshold behaviour. They can serve as a useful inside into the stochastic models as well. Even with large population, chance fluctuation may not average out, especially when these fluctuations have some spatial nature. Therefore, it may be important to take these variability under consideration. Allowance should be made for complicating feature of real infections, such as population structure and duration of disease stages - the assumed homogeneous mixing and exponential distribution of simple models are seldom appropriate. A particular challenge for the future now is to extend the epidemic model types, to allow structured population where assumptions of homogeneous mixing do not apply or partially apply. The model described in this paper is a variation on a SIR theme, with a simple and relatively tractable mathematical structure. In particular, we assume that the hosts are identical and homogeneously mixing, but not in the entire population. We divide the population of interest into subgroups/clusters, where that hosts mix homogeneously. It is straightforward to generalize this type of mixing structure to allow the distinct subgroups to have different rates of disease spreading. The groups are not isolated from one another. We allow contacts between groups that are modelled by use of transition probabilities. A particular host while traveling in space and time randomly appears in different groups - changes environment. Therefore we introduce randomness not directly through , , and , but indirectly, through the coefficients of the SIR model. These coefficients are directed by a semi-Markov process which serves as a switching process Swishchuk & Wu.⁶

The choice of a semi-Markov process is made for the purpose of generalization. In some particular situation, when the state space is finite, a Markov process can play the role of a switching process. Semi-Markovian properties have one more advantage: sojourn time must not be exponentially distributed. We note that epidemic SIR model in random media and its averaging, merging, diffusion approximation, normal deviations and stability were considered in Swishchuk & Wu.⁶ We consider an averaging principle for the endemic SIR model in a semi-Markov random media. Under stationary conditions of a semi-Markov media we show that the perturbed endemic SIR model converges to the classic endemic SIR model with averaged coefficients. Numerical toy examples and their interpretations are also presented for two-state Markov and semi-Markov chains. We also discuss two numerical examples involving real data: 1) Dengue Fever Disease (Indonesia and Malaysia (2009)) and 2) Cholera Outbreak in Zimbabwe (2008-2009). Novelty of the paper consists in studying of an endemic SIR model in semi-Markov random media and in implementations and interpretations of the results through numerical toy examples and numerical examples with real data. We note that the dengue fever model was first introduced in Derouchi et al.⁷ A SIR model for speard of dengue fever disease with simulations for South Sulawesi, Indonesia and Selanor, Malaysia, was studied in Side et al.⁸ Estimating the reprodictive numbers for the 2008-2009 cholera outbreaks in Zimbabwe was considered in Mukandavire et al.⁹ A generalized cholera model and epidemic-endemic analysis was investigated in Wang et al.¹⁰ The paper is organized as follows. Section 2 describes classic endemic SIR model. Random media is discussed in section 3. Endemic SIR model in random media is introduced in section 4. Averaged endemic SIR model is investigated in section 5. Numerical toy examples with two-state Markov and semi-Markov chains are presented in section 6. Here we also give the interpretations of the obtained theoretical and numerical results. In Section 7 we discuss two numerical examples involving real data: 1) Dengue Fever Disease (Indonesia and Malaysia (2009)) and 2) Cholera Outbreak in Zimbabwe (2008-2009). Section 8 concludes the paper and highlights some future work.

Let $S\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}I\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ and $R\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ be the number of individuals in each class of susceptible, infectives and removed, respectively. We follow the approach suggested by Hethcote⁴ taking under consideration demography and introducing $\beta$, $\mu$$\mu$, $\gamma$ as infective contact (transmission), death (mortality)/birth, and removal (or recovery) rates, respectively.

The deterministic SIR endemic model is ($N$is the number of the individuals in the population):

$\{\begin{array}{l}\frac{dS}{dt}\mathrm{<mo>=</mo>}-\beta SI+\mu \mathrm{<mo\; stretchy="false">(</mo>}N-S\mathrm{<mo\; stretchy="false">)</mo>,}S\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{S}_0\\ \\ \frac{dI}{dt}\mathrm{<mo>=</mo>}\beta SI-\gamma I-\mu I\mathrm{,}I\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{I}_0 (1)\\ \\ \frac{dR}{dt}\mathrm{<mo>=</mo>}\gamma I-\mu R\mathrm{,}R\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{R}_0\mathrm{.}\end{array}$

Sometimes it is better to work on longer time scale. Let $s\mathrm{<mo>=</mo>}S\mathrm{<mo>/</mo>}N\mathrm{,}i\mathrm{<mo>=</mo>}I\mathrm{<mo>/</mo>}N\mathrm{,}r\mathrm{<mo>=</mo>}R\mathrm{<mo>/</mo>}N\mathrm{.}$ Then the endemic SIR model in (1) becomes:

$\{\begin{array}{l}\frac{ds}{dt}\mathrm{<mo>=</mo>}-\beta si+\mu \mathrm{<mo\; stretchy="false">(</mo>1}-s\mathrm{<mo\; stretchy="false">)</mo>,}s\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{S}_0\mathrm{<mo>/</mo>}N\\ \\ \frac{di}{dt}\mathrm{<mo>=</mo>}\beta si-\gamma i-\mu i\mathrm{,}i\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{I}_0\mathrm{<mo>/</mo>}N (1)\\ \\ \frac{dr}{dt}\mathrm{<mo>=</mo>}\gamma i-\mu r\mathrm{,}r\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{R}_0\mathrm{<mo>/</mo>}N\mathrm{.}\end{array}$

In the classical endemic SIR model the various classes are uniformly mixed, that is, every pair of individuals has equal probability of coming into contact with each other and the total population. Size is constant. For many diseases with transmission taking place within some particular groups, it is logical to divide the host population into groups, where it is assumed that hosts mix homogeneously within the group. Contacts among groups are modelled by use of a transition probability matrix whose $\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{,}j\mathrm{<mo\; stretchy="false">)</mo>}$element specifies the probability that host in group $i$ will have a potential contact with host in group$j$by visiting this group. As a good example of the described above situation an epidemic on some finite amount of islands with birds populations on each can be used. Another example is diseases spreading among big cities joined with some transportation systems. If we take under consideration that one particular host may visit different groups by moving in space and time, then an important concern is how to generalize the concept of the reproduction ratio for such a heterogeneous structured population.

One of the stochastic models capable to incorporate distinct subgroups with different contact rates is the epidemic SIR model in random environment.⁶ In this model we presume that the coefficients of transmission, recovery, and mortality depend on semi-Markov process that switches the coefficient values depending on the state of the process. It looks like the system is submerged into some random media.

Let${\mathrm{<mo\; stretchy="false">(</mo>}{y}_n\mathrm{<mo\; stretchy="false">)</mo>}}_{n\in Z^{+}}$be a homogeneous Markov chain in a measurable space$\mathrm{<mo\; stretchy="false">(</mo>}Y\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>}$with transition probabilities$P\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{,}A\mathrm{<mo\; stretchy="false">)</mo>,}y\in Y\mathrm{,}A\in Y\mathrm{,}$and ergodic distribution$p\mathrm{<mo\; stretchy="false">(</mo>}A\mathrm{<mo\; stretchy="false">)</mo>,}A\in Y$;$\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}$,$\gamma \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}$, and$\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}$are non-negative, bounded measurable functions defined on$Y;{\mathrm{<mo\; stretchy="false">(</mo>}{y}_n\mathrm{<mo>;</mo>}{\theta }_n\mathrm{<mo\; stretchy="false">)</mo>}}_{n\in Z^{+}}$is a Markov renewal process in the phase space $\mathrm{<mo\; stretchy="false">(</mo>}Y\times R^{+}\mathrm{,}Y\times R^{+}\mathrm{<mo\; stretchy="false">)</mo>}$with stochastic kernel$Q\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{,}dz\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{,}dz\mathrm{<mo\; stretchy="false">)</mo>}{G}_y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}y\in Y\mathrm{,}dz\in Y,t\in R^{+}$;$\nu \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>:</mo><mo>=</mo>}\text{m}ax\mathrm{<mo>\{</mo>}n\mathrm{<mo>:</mo>}{\tau }_n\le t\mathrm{<mo>\}</mo>}$is a counting process, gis the distribution function of the sojourn times. In the discrete time case with the number of states finite, say, n, then semi-Markov stochastic kernel $Q$ has a form $Q_{ij}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{p}_{ij}{G}_i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}$where $P\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{p}_{ij}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{y}_{n+1}\mathrm{<mo>=</mo>}j\mathrm{<mo>/</mo>}{y}_n\mathrm{<mo>=</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}$is the matrix of transition probabilities, and $G_i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{\tau }_n\mathrm{<mo><</mo>}t\mathrm{<mo>/</mo>}{y}_n\mathrm{<mo>=</mo>}i\mathrm{<mo\; stretchy="false">)</mo>.}$The Markov renewal process is a convenient constructive tool to define a semi-Markov process,  ${\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}_{t\in R^{+}},y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{y}_{\nu \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}$. As$\nu \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}n,{\tau }_n\le t\mathrm{<mo><</mo>}{\tau }_{n+1},y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$also assumes constant values on the same intervals and is continuous from the right. Namely, $y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{y}_n\mathrm{,}{\tau }_n\le t\mathrm{<mo><</mo>}{\tau }_{n+1}$, and $y\mathrm{<mo\; stretchy="false">(</mo>}{\tau }_n\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{y}_n,n\ge 0.$For the semi-Markov process${\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}_{t\in R^{+}}$, the renewal time ${\theta }_n\mathrm{<mo>:</mo><mo>=</mo>}{\tau }_{n+1}-{\tau }_n$ may be naturally interpreted as the occupation time (life-time) in the state$y_n$. That explains the choice of the process: semi-Markov process differs from a Markov process by the distribution of time, for a Markov process a distribution function is exponential while for a semi-Markov process it can be any distribution function. Therefore, by choosing a semi-Markov process for a role of a switching process, we have wider possibilities for the occupation time intervals. We consider only a regular semi-Markov process, this is a process that with probability 1 has a finite number of renewals on a finite period of time. Just as in the right-continuous Markov process, the moments of jumps are regeneration points erasing the influence of the past. The only difference is that sojourn time at a point $y$has an arbitrary distribution $G_y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>},y\in Y,t\in R^{+}$which depends on the terminal state. The ergodic theorem for a semi-Markov process ${\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}_{t\in R^{+}}$states Swishchuk and Wu,⁶ that for any measured and bounded function $$f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}$$

$P\left(\frac{1}{t}{\displaystyle {\int }_0^t}f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}ds\to \widehat{f}\mathrm{,}t\to \infty \right)\mathrm{<mo>=</mo>1}$

where

$\widehat{f}\mathrm{<mo>=</mo>}\frac{1}{m}{\displaystyle {\int }_Y}m\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}dy\mathrm{<mo\; stretchy="false">)</mo>}$, $m\mathrm{<mo>=</mo>}{\displaystyle {\int }_Y}m\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}dy\mathrm{<mo\; stretchy="false">)</mo>}$, $m\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_0^{+\infty }}t{G}_y\mathrm{<mo\; stretchy="false">(</mo>}dt\mathrm{<mo\; stretchy="false">)</mo>}$

The model in semi-Markov random media/environment is defined as Swishchuk and Wu,⁶

$\{\begin{array}{l}\frac{dS}{dt}\mathrm{<mo>=</mo>}-\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}SI+\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}N-S\mathrm{<mo\; stretchy="false">)</mo>,}\\ \\ \frac{dS}{dt}\mathrm{<mo>=</mo>}-\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}SI+\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}N-S\mathrm{<mo\; stretchy="false">)</mo>,} (2)\\ \\ \frac{dR}{dt}\mathrm{<mo>=</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}I-\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}R\end{array}$

Or, in terms of new variables $s\mathrm{,}i\mathrm{,}r\mathrm{,}$ the system in (2) takes a look (1):

$\{\begin{array}{l}\frac{ds}{dt}\mathrm{<mo>=</mo>}-\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}si+\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-s\mathrm{<mo\; stretchy="false">)</mo>,}\\ \\ \frac{ds}{dt}\mathrm{<mo>=</mo>}-\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}si+\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-s\mathrm{<mo\; stretchy="false">)</mo>,} (2)\\ \\ \frac{dr}{dt}\mathrm{<mo>=</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}i-\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}r\end{array}$

It is important to note that now the coefficients$\alpha$, $\beta$, and $\gamma$ are not constant. The state of the semi-Markov process defines their value. We may say that$y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ is serving as a switching process: depending on the time this process takes its values in different states, defining corresponding to this state coefficients of disease spreading. $s\mathrm{,}i\mathrm{,}r$are random processes as well. For the finite (n-elements) set of states we have:$S\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{S}_i\mathrm{<mo\; stretchy="false">)</mo>},I\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{I}_i\mathrm{<mo\; stretchy="false">)</mo>}$and $R\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{R}_i\mathrm{<mo\; stretchy="false">)</mo>}$with$i\mathrm{<mo>=</mo>1,2,...}n$. In this case, we may consider$n$subgroups/clusters of the population under investigation. In each group we assume homogeneous mixing, but the coefficients of disease spreading may be different for different groups, taking their values as${\beta }_i$, ${\mu }_i$, and${\gamma }_i$,$i\mathrm{<mo>=</mo>1,2,...}n$.

In order to investigate the system’s (2) equilibrium we perturb this system in the following way:

$\{\begin{array}{lll}\frac{d{S}^{\epsilon }}{dt}\mathrm{<mo>=</mo>}\epsilon \mathrm{<mo\; stretchy="false">(</mo>}-\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{S}^{\epsilon }{I}^{\epsilon }+\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}N-S^{\epsilon }\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\hfill & \hfill & \hfill \\ \frac{d{I}^{\epsilon }}{dt}\mathrm{<mo>=</mo>}\epsilon \mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{S}^{\epsilon }{I}^{\epsilon }-\gamma \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{I}^{\epsilon }-\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{I}^{\epsilon }\mathrm{<mo\; stretchy="false">)</mo>}\hfill & \hfill & \hfill \\ \frac{d{R}^{\epsilon }}{dt}\mathrm{<mo>=</mo>}\epsilon \mathrm{<mo\; stretchy="false">(</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{I}^{\epsilon }-\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{R}^{\epsilon }\mathrm{<mo\; stretchy="false">)</mo>}\hfill & \hfill & \hfill \end{array}\mathrm{<mo\; stretchy="false">(</mo>3<mo\; stretchy="false">)</mo>}$

where $\epsilon$ is a small positive parameter.

Changing the time scale $t\to \frac{t}{\epsilon }$ we transform the perturbed system (3) into the system

$\{\begin{array}{l}\frac{d{S}^{\epsilon }}{dt}\mathrm{<mo>=</mo>}-\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo>/</mo>}\epsilon \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{S}^{\epsilon }{I}^{\epsilon }+\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo>/</mo>}\epsilon \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}N-S^{\epsilon }\mathrm{<mo\; stretchy="false">)</mo>}\\ \\ \frac{d{S}^{\epsilon }}{dt}\mathrm{<mo>=</mo>}-\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo>/</mo>}\epsilon \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{S}^{\epsilon }{I}^{\epsilon }+\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo>/</mo>}\epsilon \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}N-S^{\epsilon }\mathrm{<mo\; stretchy="false">)</mo>} (4)\\ \\ \frac{d{R}^{\epsilon }}{dt}\mathrm{<mo>=</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo>/</mo>}\epsilon \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{I}^{\epsilon }-\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo>/</mo>}\epsilon \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{R}^{\epsilon }\mathrm{,}\end{array}$

which can be averaged in the following way (see Swishchuk and Wu, 2003);

when $\epsilon \to 0:\mathrm{<mo\; stretchy="false">(</mo>}{S}^{\epsilon }\mathrm{,}{I}^{\epsilon }\mathrm{,}{R}^{\epsilon }\mathrm{<mo\; stretchy="false">)</mo>}\to \mathrm{<mo\; stretchy="false">(</mo>}\widehat{S}\mathrm{,}\widehat{I}\mathrm{,}\widehat{R}\mathrm{<mo\; stretchy="false">)</mo>}$, in the sence:

for any $\delta \ge 0$, ${{\displaystyle \lim }}_{\epsilon \to 0}P\mathrm{<mo>\{</mo><mo\; stretchy="false">(</mo><mo>|</mo>}{S}^{\epsilon }-\widehat{S}\mathrm{<mo>|</mo>}+\mathrm{<mo>|</mo>}{I}^{\epsilon }-\widehat{I}\mathrm{<mo>|</mo>}+\mathrm{<mo>|</mo>}{R}^{\epsilon }-\widehat{R}\mathrm{<mo>|</mo><mo\; stretchy="false">)</mo><mo>></mo>}\delta \mathrm{<mo>\}</mo><mo>=</mo>0}$

The averaged system is

$\{\begin{array}{l}\frac{d\widehat{S}}{dt}\mathrm{<mo>=</mo>}-\widehat{\beta }\widehat{S}\widehat{I}+\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}N-\widehat{S}\mathrm{<mo\; stretchy="false">)</mo>}\\ \frac{d\widehat{I}}{dt}\mathrm{<mo>=</mo>}\widehat{\beta }\widehat{S}\widehat{I}-\widehat{\gamma }\widehat{I}-\widehat{\mu }\widehat{I}\\ \frac{d\widehat{I}}{dt}\mathrm{<mo>=</mo>}\widehat{\beta }\widehat{S}\widehat{I}-\widehat{\gamma }\widehat{I}-\widehat{\mu }\widehat{I}\end{array}(5)$

In terms of new variables $s\mathrm{,}i\mathrm{,}r$ (see (1’)) the system (5) has the following form:

$\{\begin{array}{l}\frac{d\widehat{s}}{dt}\mathrm{<mo>=</mo>}-\widehat{\beta }\widehat{s}\widehat{i}+\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>1}-\widehat{s}\mathrm{<mo\; stretchy="false">)</mo>}\\ \frac{d\widehat{i}}{dt}\mathrm{<mo>=</mo>}\widehat{\beta }\widehat{s}\widehat{i}-\widehat{\gamma }\widehat{i}-\widehat{\mu }\widehat{i}\\ \frac{d\widehat{r}}{dt}\mathrm{<mo>=</mo>}\widehat{\gamma }\widehat{i}-\widehat{\mu }\widehat{r}\end{array}$

Coefficients of the averaged system can be found by using the ergodic theorem for a semi-Markov process:

$\begin{array}{l}\mu =\frac{1}{m}{\displaystyle {\int }_Y}m\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}dy\mathrm{<mo\; stretchy="false">)</mo>,}\\ \\ \beta =\frac{1}{m}{\displaystyle {\int }_Y}m\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\beta \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}dy\mathrm{<mo\; stretchy="false">)</mo>,} (6)\\ \\ \gamma =\frac{1}{m}{\displaystyle {\int }_Y}m\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}dy),\end{array}$

and

$\begin{array}{l}m={\displaystyle {\int }_Y}m\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}dy\mathrm{<mo\; stretchy="false">)</mo>,}\\ m\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}={\displaystyle {\int }_0^{+\infty }}t{G}_y\mathrm{<mo\; stretchy="false">(</mo>}dt\mathrm{<mo\; stretchy="false">)</mo>,} (7)\end{array}$

where $\pi \mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}$ is a unique invariant (stationary or ergodic) distribution. In the case when the state space of a semi-Markov process is finite,$n$-element, we have $P\mathrm{<mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{y}_{n+1}\mathrm{<mo>=</mo>}j\mathrm{<mo>/</mo>}{y}_n\mathrm{<mo>=</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>:</mo><mo>=</mo><mo\; stretchy="false">(</mo>}{p}_{ij}\mathrm{<mo>;</mo>}i\mathrm{,}j\mathrm{<mo>=</mo>1,2,...,}n\mathrm{<mo\; stretchy="false">)</mo>,}$ and the above integrals are becoming the following sums:

$\begin{array}{l}\mu =\frac{1}{m}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{m}_i\mu \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}{\pi }_i\\ \\ \beta =\frac{1}{m}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{m}_i\beta \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}{\pi }_i\mathrm{,} (8)\\ \\ \gamma =\frac{1}{m}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{m}_i\gamma \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}{\pi }_i\mathrm{,}\end{array}$

where $m_i\mathrm{<mo>=</mo>}{\displaystyle {\int }_0^{\infty }}t{G}_i\mathrm{<mo\; stretchy="false">(</mo>}dt\mathrm{<mo\; stretchy="false">)</mo>}$with a distribution function

$G_i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{\tau }_n\mathrm{<mo><</mo>}t\mathrm{<mo>/</mo>}{y}_n\mathrm{<mo>=</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}$and the stationary probabilities $\pi \mathrm{<mo>=</mo>}\left(\begin{array}{c}{\pi }_1\\ {\pi }_2\\ \mathrm{...}\\ {\pi }_n\end{array}\right)\mathrm{.}$ Here, $m\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^n}{\pi }_i{m}_i\mathrm{,}$where${\pi }_i$are ergodic distribution of $n$-state semi-Markov process.

The averaged systems $\mathrm{<mo\; stretchy="false">(</mo>5<mo\; stretchy="false">)</mo>}$ and $\mathrm{<mo\; stretchy="false">(</mo>}{5}^{\prime }\mathrm{<mo\; stretchy="false">)</mo>,}$ being deterministic, can be used for a regular analysis. For example, an endemic equilibrium solution and a basic reproductive number for the averaged system may by found by regular method. For example, $\widehat{I}\mathrm{<mo>=</mo>}\frac{\widehat{\mu }N}{\widehat{\gamma }+\widehat{\mu }}-\frac{\widehat{\mu }}{\widehat{\beta }}\mathrm{,}{\widehat{R}}_0\mathrm{<mo>=</mo>}\widehat{\beta }\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}\widehat{\gamma }+\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo>}$, etc.

In this Section, we consider a numerical toy example with two-state Markov and semi-Markov chains, and show how to find averaged data in both cases. We also explain and give insight into the results for these specific examples and show in details how the methods work.

Two-state markov chain

If$y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ is a Markov Chain with two states$y\mathrm{<mo>:</mo><mo>=</mo>}i\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>0,1<mo\; stretchy="false">)</mo>}$ and transition matrix

$P\mathrm{<mo>=</mo>}\left(\begin{array}{cc}0.7& 0.3\\ 0.4& 0.6\\ & \end{array}\right)$, the stationary probabilities now are $\overrightarrow{\pi }\mathrm{<mo>=</mo>}\left(\begin{array}{c}{\pi }_0\\ {\pi }_1\end{array}\right)\mathrm{<mo>=</mo>}\left(\begin{array}{c}0.571\\ 0.429\\ \end{array}\right)\mathrm{.}$

For the distribution function $G_i\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1}-e^{-\lambda \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}x}$(we take here exponential distribution for simplicity, but it could be taken any, e.g., gamma or Weibull (see sec. 6.3), etc.) we take

$\lambda (0)=8 and \lambda \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>10.} \mathrm{<mo\; stretchy="false">(</mo>9<mo\; stretchy="false">)</mo>}$

In this way (7),

$m\mathrm{<mo>=</mo>}{\pi }_0\mathrm{<mo>/</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}+{\pi }_1\mathrm{<mo>/</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.571<mo>/</mo>8}+\mathrm{0.429<mo>/</mo>10<mo>=</mo>0.11.} \mathrm{<mo\; stretchy="false">(</mo>10<mo\; stretchy="false">)</mo>}$

We note, that $m\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>0,1,}$ for exponential distribution $G_i\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>.}$

Suppose that our parameters $\mu \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}\gamma \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}\beta \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>0,1,}$have the following values:

$\begin{array}{l}\mu \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>60,} \mu \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>40,}\\ \gamma \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>3,} \gamma \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>4,} (11)\\ \beta \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>1.05,} \beta \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>1.1.}\end{array}$ 

Then $\widehat{\mu }\mathrm{,}\widehat{\gamma }\mathrm{,}\widehat{\beta }$ become (we use below formulas (7) and (8)):

$\begin{array}{l}\widehat{\mu }=\mathrm{0.02,}\\ \widehat{\gamma }=\mathrm{0.31,} (12)\\ \widehat{\beta }=\mathrm{1.11.}\end{array}$ 

From (12) we can find that ${\widehat{R}}_0\mathrm{<mo>=</mo>}\widehat{\beta }\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}\widehat{\gamma }+\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>3.36.}$Using these data we can find the equilibrium states $\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{s}}_e\mathrm{,} {\widehat{i}}_e\mathrm{<mo\; stretchy="false">)</mo><mo>:</mo>}$

$\begin{array}{l}{\widehat{S}}_e=\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\gamma }+\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}\widehat{\beta }\mathrm{<mo>=</mo>0.3,}\\ {\widehat{I}}_e=\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{R}}_0-\mathrm{1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>/</mo>}\widehat{\beta }\mathrm{<mo>=</mo>0.043.} (13)\end{array}$

If we take, for example, $N\mathrm{<mo>=</mo>1,000,}$then going back to our initial variables S and I we can find that equilibrium states are:

$\begin{array}{l}{\widehat{S}}_e=\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\gamma }+\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}\widehat{\beta }\times N\mathrm{<mo>=</mo>0.3}\times \mathrm{1,000<mo>=</mo>300,}\\ {\widehat{I}}_e=\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{R}}_0-\mathrm{1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>/</mo>}\widehat{\beta }\times N\mathrm{<mo>=</mo>0.043}\times \mathrm{1,000<mo>=</mo>43.} (14)\end{array}$  

Also, the value for$\widehat{I}$ is:

$\widehat{I}\mathrm{<mo>=</mo>}\widehat{\mu }\times N\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}\widehat{\mu }+\widehat{\gamma }\mathrm{<mo\; stretchy="false">)</mo>}-\widehat{\mu }\mathrm{<mo>/</mo>}\widehat{\beta }\mathrm{<mo>=</mo>61.}$ (15)

Two-state semi-markov chain

Here, we consider the case of two-state semi-Markov chain with arbitrary distribution $G_i\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>0,1,}$ for ${\tau }_n\mathrm{.}$Let us take Weibull distribution Krishnamoorthy,¹¹$G_i\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$for${\tau }_n$(see Sec. 3) with probability density function $f_i\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>:</mo><mo>=</mo>}d{G}_i\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}dx\mathrm{<mo>:</mo>}$

$f_i\left(x\right)=\{\begin{array}{l}\lambda (i)K\left(i\right){\left(\lambda \left(i\right)x\right)}^{K\left(i\right)-1}\exp \left[-{\left(\lambda \left(i\right)x\right)}^{K\left(i\right)}\right],X\le 0,\\ 0,X<0,\end{array}$(16)

where $i\mathrm{<mo>=</mo>0,1.}$ Recall that $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}$ is called shape parameter, and $\lambda \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}$-scale parameter. We note, that if we take $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1,}i\mathrm{<mo>=</mo>0,1,}$ then we have exponential distribution considered in sec. 6.2. Suppose that $\lambda \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>8}$ and $\lambda \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>10.}$We recall that the mean value for r.v. with Weibull density distribution in (16) is $\mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>1}+\mathrm{1<mo>/</mo>}K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,}$ where $\Gamma \mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$stands for Gamma distribution. We consider two cases here: i) $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>2}\mathrm{<mo\; stretchy="false">(</mo>}K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>></mo>1<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>0,1,}$and  ii) $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>2}\mathrm{<mo\; stretchy="false">(</mo>}K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo><</mo>1<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>0,1.}$The case $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1,}i\mathrm{<mo>=</mo>0,1,}$refers to exponential distribution $G_i\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$and has already been considered in sec. 6.1.

 i) Let us take $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>2}\mathrm{<mo\; stretchy="false">(</mo>}K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>></mo>1<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>0,1.}$Thus, we can calculate the following parameters, $m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}$ and $m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>:</mo>}$

$\begin{array}{l}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}=\mathrm{1<mo>/</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>1}+\mathrm{1<mo>/</mo>}K\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>1<mo>/</mo>8<mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>3<mo>/</mo>2<mo\; stretchy="false">)</mo>}\approx 0.11\\ m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}=\mathrm{1<mo>/</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>1}+\mathrm{1<mo>/</mo>}K\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>1<mo>/</mo>10<mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>3<mo>/</mo>2<mo\; stretchy="false">)</mo>}\approx \mathrm{0.088.}\end{array}$(17)

Then (see (7))

$\begin{array}{l}m={\pi }_{0}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}+{\pi }_{1}m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.571}\times 0.11+0.429\times 0.088\\ =0.06281+0.037752\\ =0.100562\approx \mathrm{0.1.}\end{array}$(18)

Suppose that the parameters $\mu \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}\beta \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}\gamma \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>.}i\mathrm{<mo>=</mo>0,1,}$are the same as in (11). Then we have for the averaged parameters:

$\begin{array}{l}\widehat{\mu }=\frac{1}{m}\mathrm{<mo\; stretchy="false">[</mo>}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\pi }_0+m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}{\pi }_1\mathrm{<mo\; stretchy="false">]</mo>}\\ \\ =\frac{1}{0.1}\mathrm{<mo\; stretchy="false">[</mo>0.11}\times \mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>60<mo\; stretchy="false">)</mo>}\times 0.571+0.088\times \mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>40<mo\; stretchy="false">)</mo>}\times \mathrm{0.429<mo\; stretchy="false">]</mo>}\\ \\ =\mathrm{10<mo\; stretchy="false">[</mo>0.0010468}+\mathrm{0.0009438<mo\; stretchy="false">]</mo><mo>=</mo>0.019906}\approx 0.02\\ \widehat{\gamma }=\frac{1}{m}\mathrm{<mo\; stretchy="false">[</mo>}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\pi }_0+m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}{\pi }_1\mathrm{<mo\; stretchy="false">]</mo>}\\ \\ =\frac{1}{0.1}\mathrm{<mo\; stretchy="false">[</mo>0.11}\times \mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>3<mo\; stretchy="false">)</mo>}\times 0.571+0.088\times \mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>4<mo\; stretchy="false">)</mo>}\times \mathrm{0.429<mo\; stretchy="false">]</mo>}\\ \\ =\mathrm{10<mo\; stretchy="false">[</mo>0.0209366}+\mathrm{0.009438<mo\; stretchy="false">]</mo><mo>=</mo>0.0303746}\approx 0.03\\ \widehat{\beta }=\frac{1}{m}\mathrm{<mo\; stretchy="false">[</mo>}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\beta \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\pi }_0+m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}\beta \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}{\pi }_1\mathrm{<mo\; stretchy="false">]</mo>}\\ \\ =\frac{1}{0.1}\mathrm{<mo\; stretchy="false">[</mo>0.11}\times \mathrm{<mo\; stretchy="false">(</mo>1.05<mo\; stretchy="false">)</mo>}\times 0.571+0.088\times \mathrm{<mo\; stretchy="false">(</mo>1.1<mo\; stretchy="false">)</mo>}\times \mathrm{0.429<mo\; stretchy="false">]</mo>}\\ \\ =\mathrm{10<mo\; stretchy="false">[</mo>0.0659505}+\mathrm{0.0415272<mo\; stretchy="false">]</mo><mo>=</mo>0.1074777}\approx \mathrm{0.11,}\end{array}$ (19)

where$m$was calculated in (18). From here we can find

${\widehat{R}}_0\mathrm{<mo>=</mo>}\widehat{\beta }\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}\widehat{\gamma }+\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>0.11<mo>/</mo><mo\; stretchy="false">(</mo>0.03}+\mathrm{0.02<mo\; stretchy="false">)</mo><mo>=</mo>0.11<mo>/</mo>0.05<mo>=</mo>2.2}$ (20)

Also, the equilibrium states are:

$\begin{array}{l}{\widehat{s}}_e=\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\gamma }+\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}\widehat{\beta }\\ =\mathrm{0.05<mo>/</mo>0.11<mo>=</mo>0.004545}\\ =\approx 0.0045\\ {\widehat{i}}_e=\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{R}}_0-I\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>/</mo>}\widehat{\beta }\\ =\mathrm{0.02<mo\; stretchy="false">(</mo>2.2}-\mathrm{1<mo\; stretchy="false">)</mo><mo>/</mo>0.11<mo>=</mo>0.02}\times 1.2\times \mathrm{0.11<mo>=</mo>0.00264}\\ =\approx \mathrm{0.003.}\end{array}$ (21)

If we take again$N\mathrm{<mo>=</mo>1,000,}$then from (21) we have:

$\begin{array}{l}{\widehat{S}}_e={\widehat{s}}_e\times \mathrm{1,000<mo>=</mo>4.5}\\ {\widehat{I}}_e={\widehat{i}}_e\times \mathrm{1,000<mo>=</mo>0.003}\times \mathrm{1,000<mo>=</mo>3.}\end{array}$(22)

Also, the value for$\widehat{I}$is:

$\begin{array}{l}\widehat{I}=\widehat{\mu }\times N\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}\widehat{\mu }+\widehat{\gamma }\mathrm{<mo\; stretchy="false">)</mo>}-\widehat{\mu }\mathrm{<mo>/</mo>}\widehat{\beta }\\ =0.02\times \mathrm{1,000<mo>/</mo><mo\; stretchy="false">(</mo>0.02}+\mathrm{0.03<mo\; stretchy="false">)</mo>}-\mathrm{0.02<mo>/</mo>0.11}\\ =0.02\times \mathrm{1,000<mo>/</mo><mo\; stretchy="false">(</mo>0.02}+\mathrm{0.03<mo\; stretchy="false">)</mo>}-\mathrm{0.02<mo>/</mo>0.11}\end{array}$ (23)

 ii) Let us take $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>2}\mathrm{<mo\; stretchy="false">(</mo>}K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo><</mo>1<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>0,1.}$

Thus, we can again calculate the following parameters, $m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}$and $m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}$

$\begin{array}{l}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}=\mathrm{1<mo>/</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>1}+\mathrm{1<mo>/</mo>}K\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>1<mo>/</mo>8<mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>3<mo\; stretchy="false">)</mo>}\approx 0.25\\ m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}=\mathrm{1<mo>/</mo>}\lambda \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>1}+\mathrm{1<mo>/</mo>}K\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo><mo\; stretchy="false">(</mo>1<mo>/</mo>10<mo\; stretchy="false">)</mo>}\Gamma \mathrm{<mo\; stretchy="false">(</mo>3<mo\; stretchy="false">)</mo>}\approx \mathrm{0.2.}\end{array}$ (24) Then (see (7))

$\begin{array}{l}m={\pi }_{0}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}+{\pi }_{1}m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.571}\times 0.25+0.429\times 0.2\\ =0.14275+\mathrm{0.0858<mo>=</mo>0.22855}\approx 0.23\end{array}$ (25)

Suppose that the parameters$\mu \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}\beta \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}\gamma \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>.}i\mathrm{<mo>=</mo>0,1,}$are the same as in (11). Then we have for the averaged parameters (see (24)):

$\begin{array}{l}\widehat{\mu }=\frac{1}{m}\mathrm{<mo\; stretchy="false">[</mo>}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\pi }_0+m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}{\pi }_1\mathrm{<mo\; stretchy="false">]</mo>}\\ \\ =\frac{1}{0.23}\mathrm{<mo\; stretchy="false">[</mo>0.25}\times \mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>60<mo\; stretchy="false">)</mo>}\times 0.571+0.2\times \mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>40<mo\; stretchy="false">)</mo>}\times \mathrm{0.429<mo\; stretchy="false">]</mo>}\\ \\ \widehat{\gamma }=\frac{1}{m}\mathrm{<mo\; stretchy="false">[</mo>}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\pi }_0+m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}\gamma \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}{\pi }_1\mathrm{<mo\; stretchy="false">]</mo>}\\ \\ =\frac{1}{0.23}\mathrm{<mo\; stretchy="false">[</mo>0.25}\times \mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>3<mo\; stretchy="false">)</mo>}\times 0.571+0.2\times \mathrm{<mo\; stretchy="false">(</mo>1<mo>/</mo>4<mo\; stretchy="false">)</mo>}\times \mathrm{0.429<mo\; stretchy="false">]</mo>}\\ \\ =\mathrm{4.35<mo\; stretchy="false">[</mo>0.08}+\mathrm{0.02<mo\; stretchy="false">]</mo><mo>=</mo>0.435}\approx 0.44\\ \widehat{\beta }=\frac{1}{m}\mathrm{<mo\; stretchy="false">[</mo>}m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\beta \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\pi }_0+m\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}\beta \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}{\pi }_1\mathrm{<mo\; stretchy="false">]</mo>}\\ \\ =\frac{1}{0.23}\mathrm{<mo\; stretchy="false">[</mo>0.25}\times \mathrm{<mo\; stretchy="false">(</mo>1.05<mo\; stretchy="false">)</mo>}\times 0.571+0.2\times \mathrm{<mo\; stretchy="false">(</mo>1.1<mo\; stretchy="false">)</mo>}\times \mathrm{0.429<mo\; stretchy="false">]</mo>}\\ \\ =\mathrm{4.35<mo\; stretchy="false">[</mo>0.15}+\mathrm{0.09<mo\; stretchy="false">]</mo>}\approx \mathrm{1.04,}\end{array}$ (26)

where$m$was calculated in (24). From here we can find

${\widehat{R}}_0\mathrm{<mo>=</mo>}\widehat{\beta }\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}\widehat{\gamma }+\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1.04<mo>/</mo><mo\; stretchy="false">(</mo>0.44}+\mathrm{0.02<mo\; stretchy="false">)</mo><mo>=</mo>1.04<mo>/</mo>0.4}\approx 2.26$(27)

Also, the equilibrium states are:

$\begin{array}{l}{\widehat{s}}_e=\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\gamma }+\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}\widehat{\beta }\\ =\mathrm{0.46<mo>/</mo>1.04}\\ =0.44\\ {\widehat{i}}_e=\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\mu }\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{R}}_0-I\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>/</mo>}\widehat{\beta }\\ =\mathrm{0.02<mo\; stretchy="false">(</mo>2.26}-\mathrm{1<mo\; stretchy="false">)</mo><mo>/</mo>1.04}\\ =0.02\end{array}$(28)

If we take again $N\mathrm{<mo>=</mo>1,000,}$ then from (27) we have:

$\begin{array}{l}{\widehat{S}}_e={\widehat{\widehat{s}}}_e\times \mathrm{1,000<mo>=</mo>0.44}\times \mathrm{1,000<mo>=</mo>440}\\ {\widehat{I}}_e={\widehat{i}}_e\times \mathrm{1,000<mo>=</mo>0.0.02}\times \mathrm{1,000<mo>=</mo>20.}\end{array}$ (29)

Also, the value for$\widehat{I}$is:

$\begin{array}{l}\widehat{I}=\widehat{\mu }\times N\mathrm{<mo>/</mo><mo\; stretchy="false">(</mo>}\widehat{\mu }+\widehat{\gamma }\mathrm{<mo\; stretchy="false">)</mo>}-\widehat{\mu }\mathrm{<mo>/</mo>}\widehat{\beta }\\ =0.02\times \mathrm{1,000<mo>/</mo><mo\; stretchy="false">(</mo>0.46<mo\; stretchy="false">)</mo>}-\mathrm{0.02<mo>/</mo>1.04}\\ =\mathrm{20<mo>/</mo>0.46}-\mathrm{0.019<mo>=</mo>43.482}-\mathrm{0.019<mo>=</mo>43.4692}\\ =43.\end{array}$

Interpretation of the numerical examples

In our case, we have an infection that is endemic in a community when transmission persists. Including two-state Markov chain into coefficients $\mu \mathrm{,}\gamma \mathrm{,}\beta$ means that the endemic develops with respect to two modes: one mode $0$ with coefficients $\mu \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>,}\gamma \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>,}\beta \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}$and another mode 1 with coefficients $\mu \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>,}\gamma \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>,}\beta \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}$(see (11)). Probability to stay in mode 0 is 0.7, probability to stay in mode 0 is 0.7 probability to switch from mode 0 to mode 1 is 0.3 and probability to switch from mode 1 to mode 0 is 0.4 (see entries of matrix P).

Two-state markov chain case: sec. 6.1

Time to stay in mode 0 or mode 1 distributed exponentially with parameter $\lambda \mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>,}i\mathrm{<mo>=</mo>0,1,}$respectively. Initially, the contact numbers were $R_0^0\mathrm{<mo>=</mo>}\beta \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>/</mo><mo\; stretchy="false">(</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}+\gamma \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo>3}$for mode 0,and $R_0^1\mathrm{<mo>=</mo>}\beta \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>/</mo><mo\; stretchy="false">(</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}+\gamma \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo>4}$for mode 1 After averaging it is ${\widehat{R}}_0\mathrm{<mo>=</mo>3.36}$(see sec. 6.1). The solution to the endemic averaged SIR endemic model (5’) eventually settles down to a steady state. We defined this steady state $\mathrm{<mo\; stretchy="false">(</mo>}{s}_e\mathrm{,}{i}_e\mathrm{<mo\; stretchy="false">)</mo>}$in (13) by solving the equations ${\widehat{s}}^{\prime }\mathrm{<mo>=</mo>0}$and ${\widehat{i}}^{\prime }\mathrm{<mo>=</mo>0.}$It means that on the long time interval the number ${\widehat{S}}_e$of susceptible is ${\widehat{S}}_e\mathrm{<mo>=</mo>300,}$and the number ${\widehat{I}}_e$of infectious is ${\widehat{I}}_e\mathrm{<mo>=</mo>43}$(see (14)). If we compare the latter number ${\widehat{I}}_e\mathrm{<mo>=</mo>43}$(steady state for $\widehat{I}$) with $\widehat{I}\mathrm{<mo>=</mo>61}$then we can see that the number of infectious is decreasing on the long time interval.

Two-state semi-markov chain case: sec. 6.2

The parameter$K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo>}$in Weibull distribution describes the  failure rate for the disease: if $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo><</mo>1,}$then the number of infectious of the disease decreases over time, if $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1,}$then the number of infectious of the disease is constant over time (exponential distribution), and if $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>></mo>1,}$then the number of infectious of the disease increases over time. In our cases: if  i) $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>2<mo>></mo>1,}i\mathrm{<mo>=</mo>0,1,}$then it means that the number of infectious of the disease in this case increases. And it was really the case: the number of infectious $\widehat{I}$increased to 400 (see (23)), but after the long-time period decreased and stabilized to 3 (see (22)). Compare with the Markov case, where $\widehat{I}\mathrm{<mo>=</mo>61,}$the semi-Markov case gave us much bigger number of infectious, namely $\widehat{I}\mathrm{<mo>=</mo>400.}$if ii)$K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>2<mo><</mo>1,}i\mathrm{<mo>=</mo>0,1,}$then it means that the number of infectious of the disease in this case decreases. And it was really the case: the number of infectious$\widehat{I}$ was 43 (see (30)), but after the long-time period decreased and stabilized to 20 (see (29)). Compare with the Markov case, where $\widehat{I}\mathrm{<mo>=</mo>61,}$the semi-Markov case gave us much lower number of infectious, namely $\widehat{I}\mathrm{<mo>=</mo>43.}$Thus, the result crucially depends on a distribution of time the chain spent in a state. In our case, we compared two distributions, exponential and Weibull. Similar numerical examples can be prepared for other distributions $G_i\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>,}$such as Gamma or Beta, etc. Krishnamoorthy.^11,12

In this section we discuss two numerical examples involving real data: 1) Dengue Fever Disease (Indonesia and Malaysia (2009)) and 2) Cholera Outbreak in Zimbabwe (2008-2009). We show how to construct two-state Markov and semi-Markov chains and to obtain some estimations associated with these real data. In the first case we use two countries, Indonesia and Malaysia, and in the second case we use one country, Zimbabwe, but with many (in fact, 10) regions, and take two regions to get the two state Markov or semi-Markov chains. These real data sets were borrowed from respectively.^8,9 Crucial problems in the real data examples are: 1) determine the matrix P of transition probabilities between different countries or different regions for one country, which are the states of our Markov or semi-Markov chains; 2) determine whether we have Markov or semi-Markov chain.

Spread of dengue fever disease (south sulawesi, indonesia, and selangor, malaysia (2009))

As long as the dengue fever disease was spread in two countries, Indonesia and Malaysia, we can consider these two countries as two states of our Markov or semi-Markov chains. Let us take state ’0’ for Indonesia and state ’1’ for Malaysia. As long as we do not have any data for migration between those two countries, which we need to create our matrix P we will use our matrix from sec. 6.1:$P\mathrm{<mo>=</mo>}\left(\begin{array}{cc}0.7& 0.3\\ 0.4& 0.6\end{array}\right)$,

with the stationary probabilities $\overrightarrow{\pi }\mathrm{<mo>=</mo>}\left(\begin{array}{c}{\pi }_0\\ {\pi }_1\end{array}\right)\mathrm{<mo>=</mo>}\left(\begin{array}{c}0.571\\ 0.429\\ \end{array}\right)\mathrm{.}$

If we had the migration information then the matrix $P$could be calculated easily using law of large numbers. Probabilities in matrix$P$can be interpreted as follows: 0.7 means probability to stay in Indonesia, 0.3-probability to move from Indonesia to Malaysia, 0.4 -probability to move from Malaysia to Indonesia,0.6 and -probability to stay in Malaysia. Meaning of stationary probabilities: 0.571 -probability to stay in Indonesia on the long time interval,0.429 -probability to stay in Malaysia on the long time interval. Further, again, as long as we do not know the migration information we cannot judge about the intensity of the migration, where it follows exponential distribution or any non-exponential distribution. If we knew this migration information then we could calculate/calibrate the intensity of migration and make a decision about the parameters of exponential distribution or non-exponential. Now, we will take the data from Side et al.⁸ for Indonesia (see page 101): $\mu \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0.000046,}\gamma \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0.328833,}\beta \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0.75<mo>;</mo>}$for Malaysia (see page 102): $\mu \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.0045,}\gamma \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.15,}\beta \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.75.}$We consider, again, two cases: 1) migration intensity follows exponential distribution and ii) migration intensity follows non-exponential distribution. We take again Weibull distribution in case 2) because we have already shown how to use it in sec. 6.2.

Exponential case

We suppose that$\lambda \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>10}$ and $\lambda \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>10.}$We note, that$m\mathrm{<mo>=</mo>0.1}$ then. Then our averaged parameters have the following values (see sec. 6.1 for calculations formulas):$\widehat{\mu }\mathrm{<mo>=</mo>0.00183,}\widehat{\gamma }\mathrm{<mo>=</mo>0.258,}\widehat{\beta }\mathrm{<mo>=</mo>0.75.}$The the averaged reproductive number is ${\widehat{R}}_0\mathrm{<mo>=</mo>2.8<mo>></mo>1,}$and we have the case of endemic situation. Thus, if the intensity of migration is high then we can get the endemic situation.To avoid it, we have to restrict migration between those two countries during epidemic of dengue fever disease.

Non-exponential case: weibull distribution

For the Weibull distribution (see sec. 6.2) of intensity of migration we consider again two cases:  i)$K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>2}$and  ii) $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>2,}i\mathrm{<mo>=</mo>0,1.}$Using the real data from Side  et al. (2013): for Indonesia (see page 101): $\mu \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0.000046,}\gamma \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0.328833,}\beta \mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0.75<mo>;</mo>}$ for Malaysia (see page 102): $\mu \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.0045,}\gamma \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.15,}\beta \mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo><mo>=</mo>0.75,}$we have (see sec. 6.2 for calculation formulas): i) $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>2,}i\mathrm{<mo>=</mo>0,1.}$ We note that $m\mathrm{<mo>=</mo>0.88}$ here. Then, the averaged numbers are $\widehat{\mu }\mathrm{<mo>=</mo>0.0001827,}\widehat{\gamma }\mathrm{<mo>=</mo>0.0255,}\widehat{\beta }\mathrm{<mo>=</mo>0.075.}$Then, the averaged reproductive number ${\widehat{R}}_0\mathrm{<mo>=</mo>3<mo>></mo>1,}$and we have endemic situation again. ii) $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>2,}i\mathrm{<mo>=</mo>0,1.}$We note that $m-0.2$ here. Then, the averaged numbers are$\widehat{\mu }\mathrm{<mo>=</mo>0.00036552,}\widehat{\gamma }\mathrm{<mo>=</mo>0.258,}\widehat{\beta }\mathrm{<mo>=</mo>0.75.}$ Then, the averaged reproductive number${\widehat{R}}_0\mathrm{<mo>=</mo>2.1<mo>></mo>1,}$and we have endemic situation again, even when $K\mathrm{<mo\; stretchy="false">(</mo>}i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1<mo>/</mo>2<mo><</mo>1.}$In all three cases we have endemic situation. To avoid it, we must restrict migration between those two countries.

Cholera outbreak in zimbabwe (2008-2009)

In this case of cholera outbreak in Zimbabwe, we use the real data from Mukandavire et al.⁹ We have 10 provinces in Zimbabwe with different population size, total infected sizes, attack rates per $\mathrm{10,000,}$ total deaths ( page 8769) and estimates for the basic reproductive numbers $R_0$for all 10 provinces (8th row, page 8769). In all 10 cases the reproductive numbers are greater than 1, meaning that endemicity is possible. Now, if we suppose that there is migration between 10 provinces then the situation does not change. Again, we do not have any migration information between 10 provinces, thus we can use our matrix $P$for transition probabilities between any of two provinces and stationary probabilities $\overrightarrow{\pi }$ from previous sections again. Suppose that we take 2 provinces, Mashonaland East and Balawayo, with the smallest reproductive numers,$1.11$and 1.36 respectively. Then, on the long time interval, the averaged reproductive number is ${\widehat{R}}_0\mathrm{<mo>=</mo>1.21<mo>></mo>1,}$ meaning endemicity again. The situation is similar for the rest of the provinces. To reduce the reproductive number below 1, we must restrict the migration between provinces thast have cholera outbreak.

Remark

We considered only the case of two-state Markov or semi-Markov chains. However, we could consider the case of three- or more states Markov or semi-Markov chains, and pursue with analogical calculations to get the averaged results.

In this paper, we considered a random media as a semi-Markov process because in some cases the distribution functions of being in states (for models, for example, for birds migrating between some islands, or people migrating between several cities, etc.) are not exponentially distributed.⁶ From the other side, it is more general mathematical model than Markov model. The main result of the paper is the averaging principle for endemic SIR model in semi-Markov random media. Under stationary conditions of semi-Markov media we shown that the perturbed endemic SIR model converges to the classic SIR model with averaged coefficients. We also considered numerical examples with two-state Markov and semi-Markov chains and gave interpretation of the obtained results. We also discussed two numerical examples involving real data: 1) Dengue Fever Disease (Indonesia and Malaysia (2009)) and 2) Cholera Outbreak in Zimbabwe (2008-2009). We have shown how to construct two-state Markov and semi-Markov chains and to obtain some estimations associated with these real data. In the first case we used two countries, Indonesia and Malaysia, and in the second case we use one country, Zimbabwe, but with two regions, to get the two state Markov or semi-Markov chains. This research paper is just a first step in the investigation of endemic SIR model in random media. The future work will be devoted to the merging, diffusion approximation, normal deviations and stability of endemic SIR models in semi-Markov random media.

None.

The authors declare that there is no conflict of interests regarding the publication of this paper.

1. Nasel I. The threshold concept in stochastic epidemic and endemic models. In Epidemic Models: Their Structure and Relation to data. USA: Cambridge University Press, 1985;71–83.
2. Nasel I. Stochastic models of some endemic infections. Math Biosci. 2002;179(1):1–19.
3. Ball F, Lyne O. Stochastic multi-type SIR epidemics among a population partitioned into households. Advances in Applied Probability. 2001;33(1):99–123.
4. Hethcote HW. The mathematics of infectious diseases. SIAM review 2000;42(4):599–653.
5. Allen L, Burguin A. Comparison of deterministic and stochastic SIS and SIR models in discrete time. Math Biosci. 2000;163(1):1–23.
6. Swishchuk A, Wu J. Evolution of Biological Systems in Random Media: Limit Theorems and Stability. USA: Kluwer Academic Publishers, 2003.
7. Derouchi M, Boutayeb A, Twizell EH, et al. A model of dengue fever. BioMedical Engineering Online. 2003;2(4).
8. Side S, Noorani S. A SIR model for spread of dengue fever disease (Simulation for South Sulawesi, Indonesia and Selangor, Malaisia). World J Model & Simul. 2013;9(2):96–105.
9. Mukandavire Z, Liao S, Wang J, et al. Estimating the reproductive numbers for the 2008-2009 cholera outbreaks in Zimbabwe. PNAS 2011;108(21):8767–8772.
10. Wang J, Liao S. A generalized cholera model and epidemic-endemic analysis. J Biol Dyn. 2012;62:568–589.
11. Krishnamoorthy K. Handbook of Statistical Distributions with Applications. USA: Chapman & Hall/CRC Press; 2006.
12. Isham V, Medley G. Models for Infectious Human diseases: their Structure and Relation to Data. USA: Cambridge University Press; 1996.