ABSTRACT

A simulation model based on nonlinear ordinary differential equations to interpret the transmission dynamics of Zika Virus (ZIKV), is formulated and analyzed, integrating the asymptomatic human population and coupled to the Aedes aegypti dynamics, the epidemic threshold Basic Reproduction Number R₀ is determined, as the spectral radius of Next-Generation Matrix and the system is simulated with MAPLE computing program taking the parameter values from literature.

Keywords:

Simulation Models, Zika Virus (ZIKV), Basic Reproductive Number R₀, Spectral Radius, Asymptomatic, Next-Generation Matrix

1. Introduction

The Zika virus (ZIKV) is an arthropod-borne virus belonging to flavivirus family, is transmitted by mosquitoes of the Aedes genus [1] . It was isolated from Aedes africanus in 1948 [2] but its first appearance dates back to 1947, presenting in a monkey from the forests of Zika in Uganda, due to this its name [3] . The first manifestations of this virus in humans occurred in 1952 in East Africa. Between 1977 and 1978 there was an outbreak in Malaysia and Indonesia attributed to a strong rain season. Until 2007 on the Yap Island in Micronesia new cases were detected [4] [5] [6] .

ZIKV can be transmitted both sexually and vectorially. Sexually transmitted cases have been reported [1] [7] [8] . Contacting with fluids such as human saliva, blood, semen and urine are among the nonvectorial forms for transmission of the virus [7] . The ZIKV transmission through sexual contact has been suggested by Foy et al., who describes two patients infected in the southeast Senegal in 2008 [5] .

Until recently, the disease was characterized as clinically mild. The primary symptoms last approximately a week. Around 75% of infections belong to asymptomatic population [9] [10] , there could be a great deal of under reporting of cases since up to 80% of cases are asymptomatic [11] [12] , individuals who are not infectious (asymptomatic period) [13] and rarely there is a need of hospitalization [14] . While the disease can be asymptomatic or cause minor fever and a rash, the virus has also been linked to Guillian-Barré syndrome and microcephaly in infants born to mothers infected during pregnancy [11] . Fever, joint pains, headaches and conjunctivitis are the most relevant symptoms [15] .

In ZIKV disease, it’s important to consider the asymptomatic people due to the high probability of disseminating the virus among susceptible population in unnoticed way. Some of the studies that have been focused to the ZIKV, such as Fréour et al. (2016), who studied the sexual transmission of this disease, probably among heterosexuals, in a totally asymptomatic couple [16] . Ascertainment bias cannot be ignored for this epidemiological scourge with substantial fraction of asymptomatic and mild infections [17] . A mathematical analysis suggests that ZIKV may exhibit similar dynamics to dengue virus in island populations, with transmission characterized by large, sporadic outbreaks with a high proportion of asymptomatic or unreported cases [18] . Pacheco et al. (2016), presented a study on the general context of ZIKV in Colombia, demonstrating the relationship between the virus and microcephaly. At the same time, he gave a report about asymptomatic people to this illness [19] . When the reproduction number varies between 1.9 and 2.8 based on estimates for Antioquia-Colombia, the attack rate varies between 4% and 26% for a low (below 10%) effectiveness of interventions in blunting the ZIKV transmission from symptomatic cases to mosquitoes [20] .

On the other hand, Manrique et al. (2016) presented a mathematical model where the relationship in the transmission of the virus in two ways was explored, both sexual and vectorial, with special emphasis on the consequences of the nonvectorial transmission in the population [21] . Agusto et al. (2017) analyze a new system of ordinary differential equations, which incorporates human vertical transmission of ZIKV, the birth of babies with microcephaly and asymptomatically infected individuals. Their results show that asymptomatic individuals amplify the disease burden in the community and the model was locally and globally asymptotically stable when the R₀ was less than unity [22] . The same authors [23] , set out another model incorporating both vector and sexual transmission routes, finding when there is no disease induced mortality the disease- free equilibrium is locally and globally asymptotically stable whenever the associated reproduction number is less than unity and unstable otherwise. Padmanabhan et al. (2017), consider two classes of human infectious, asymptomatic and symptomatic, which are assumed to be equally infectious and of similar duration (period of infectiousness) [24] . Onuorah et al. (2016) suggest a deterministic mathematical model of ZIKV where obtained the R₀ using the next generation method [25] , indicating to be most sensitive to the biting and mortality rate of mosquitoes while sexual transmission increases the risk of infection and epidemic size and prolongs the outbreak [26] .

This article suggests and analyzes a simulation model based on nonlinear ordinary differential equations, which interprets the ZIKV transmission, integrating the asymptomatic population and coupled to the A. aegypti dynamics. Subsequently, in its analysis the epidemic threshold basic reproduction number R₀ by the theory of Next-Generation Matrix, is determined [25] . Finally, simulations of the differential equation system implementing MAPLE software [27] are made.

2. The Model

A model for the transmission of the ZIKV with symptomatic and asymptomatic population is formulated (see Figure 1), coupling the carrier and noncarrier adult mosquito population dynamics, considering both infected human populations may transmit ZIKV to the noncarrier mosquitoes. Constant fluxes from susceptible people and from noncarrier mosquitoes, constant rate of human mortality, natural death rate of mosquitoes by environmental conditions and incidence type standard, are assumed. The host-vector infectious process using

nonlinear ordinary differential equations is made, and the epidemic threshold basic reproduction number (i.e. the average number of new cases of infected people during infectious period in the susceptible population) by operated matricial method also known as the Next-Generation Matrix and its respective spectral radius, is determined.

Population variables are time functions, considering $x_{\text{1}}\left(t\right)=x_{\text{1}}$ , $x_{\text{2}}\left(t\right)=x_{\text{2}}$ , $x_{\text{3}}\left(t\right)=x_{\text{3}}$ , $x_{\text{4}}\left(t\right)=x_{\text{4}}$ , $y_{\text{1}}\left(t\right)=y_{\text{1}}$ y $y_{\text{2}}\left(t\right)=y_{\text{2}}$ .

Variables, initial conditions and parameters of the dynamics are shown in Table 1 and Table 2. Thus, the dynamic system that interprets the infectious process is as follows:

$\frac{\text{d}{x}_1}{\text{d}t}=\delta -\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1-\mu x_1$ (1)

$\frac{\text{d}{x}_2}{\text{d}t}=f\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1-\left(\theta +\mu \right)x_2$ (2)

$\frac{\text{d}{x}_3}{\text{d}t}=\left(1-f\right)\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1-\left(\omega +\mu \right)x_3$ (3)

$\frac{\text{d}{x}_4}{\text{d}t}=\theta x_2+\omega x_3-\mu x_4$ (4)

$\frac{\text{d}{y}_1}{\text{d}t}=\eta -\sigma \frac{y_1}{{{\displaystyle \sum }}_{i=1}^4{x}_i}{\displaystyle \sum }_{k=2}^3{x}_k-ϵy_1$ (5)

$\frac{\text{d}{y}_2}{\text{d}t}=\sigma \frac{y_1}{{{\displaystyle \sum }}_{i=1}^4{x}_i}{\displaystyle \sum }_{k=2}^3{x}_k-ϵy_2$ (6)

The equation system (1)-(6) belongs to a dynamic system where each differential equation states the variability of a population in the time. Initial conditions are: $x_{\text{1}}\left(0\right)=x_{\text{1}0}$ , $x_{\text{2}}\left(0\right)=x_{\text{2}0}$ , $x_{\text{3}}\left(0\right)=x_{\text{3}0}$ , $x_{\text{4}}\left(0\right)=x_{\text{4}0}$ , $y_{\text{1}}\left(0\right)=y_{\text{1}0}$ y $y_{\text{2}}\left(0\right)=y_{\text{2}0}$ .

We define the epidemiological meaning region, where the solution curves or trajectories of each population are positive, as:

$\text{Ω}=\left\{(x_1,x_2,x_3,x_4,y_1,y_2)\in R_{+}^6:0<{\displaystyle \sum }_{i=1}^4{x}_i\le \frac{\delta }{\mu },0\le y_1+y_2\le \frac{\eta }{ϵ}\right\}$

3. Basic Reproduction Number, R₀

To calculate the basic reproduction number R₀, as the spectral radius of the Next-Generation Matrix ρ(G) [28] [29] [30] , we rewrite the system (1) to (6), according to the equations of infectious process:

$\frac{\text{d}{x}_2}{\text{d}t}=f\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1-\left(\theta +\mu \right)x_2$

$\frac{\text{d}{x}_3}{\text{d}t}=\left(1-f\right)\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1-\left(\omega +\mu \right)x_3$

$\frac{\text{d}{y}_2}{\text{d}t}=\sigma \frac{y_1}{{{\displaystyle \sum }}_{i=1}^4{x}_i}{\displaystyle \sum }_{k=2}^3{x}_k-ϵy_2$

$\frac{\text{d}{x}_1}{\text{d}t}=\delta -\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1-\mu x_1$

$\frac{\text{d}{x}_4}{\text{d}t}=\theta x_2+\omega x_3-\mu x_4$

$\frac{\text{d}{y}_1}{\text{d}t}=\eta -\sigma \frac{y_1}{{{\displaystyle \sum }}_{i=1}^4{x}_i}{\displaystyle \sum }_{k=2}^3{x}_k-ϵy_1$

This system has the $\dot{X}=\Gamma \left(X\right)-\Phi \left(X\right)$ form, where,

$\text{Γ}\left(X\right)=\left(\begin{array}{c}f\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1\\ \left(1-f\right)\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1\\ \sigma \frac{{{\displaystyle \sum }}_{k=2}^3{x}_k}{{{\displaystyle \sum }}_{i=1}^4{x}_i}{y}_1\\ 0\\ 0\\ 0\end{array}\right)$ y $\text{Φ}\left(X\right)=\left(\begin{array}{c}\left(\theta +\mu \right)x_2\\ \left(\omega +\mu \right)x_3\\ ϵy_2\\ -\delta +\beta \frac{y_2}{{{\displaystyle \sum }}_{j=1}^2{y}_j}{x}_1+\mu x_1\\ -\theta x_2-\omega x_3+\mu x_4\\ -\eta +\sigma \frac{{{\displaystyle \sum }}_{k=2}^3{x}_k}{{{\displaystyle \sum }}_{i=1}^4{x}_i}{y}_1+ϵy_1\end{array}\right)$ (7)

With F in the equilibrium point $E_0=\left(\frac{\delta }{\mu },0,0,0,\frac{\eta }{ϵ},0\right)$ ,

$F=\left(\frac{\partial \Gamma \left(X\right)}{\partial X}\right)\left(E_0\right)=\left(\begin{array}{ccc}0& 0& \frac{f\beta {\widehat{x}}_1{\widehat{y}}_1}{{\left({{\displaystyle \sum }}_{j=1}^2{y}_j\right)}^2}\\ 0& 0& \frac{\left(1-f\right)\beta {\widehat{x}}_1{\widehat{y}}_1}{{\left({{\displaystyle \sum }}_{j=1}^2{y}_j\right)}^2}\\ \frac{\sigma {\widehat{y}}_1\left({\widehat{x}}_1+{\widehat{x}}_2\right)}{{\left({{\displaystyle \sum }}_{i=1}^4{x}_i\right)}^2}& \frac{\sigma {\widehat{y}}_1\left({\widehat{x}}_1+{\widehat{x}}_2\right)}{{\left({{\displaystyle \sum }}_{i=1}^4{x}_i\right)}^2}& 0\end{array}\right)$ (8)

$=\left(\begin{array}{ccc}0& 0& \frac{f\beta {\widehat{x}}_1}{{\widehat{y}}_1}\\ 0& 0& \frac{p{\widehat{x}}_1}{{\widehat{y}}_1}\\ \frac{\sigma {\widehat{y}}_1}{{\widehat{x}}_1}& \frac{\sigma {\widehat{y}}_1}{{\widehat{x}}_1}& 0\end{array}\right)=\left(\begin{array}{ccc}0& 0& \frac{f\beta \delta ϵ}{\mu \eta }\\ 0& 0& \frac{p\beta \delta ϵ}{\mu \eta ϵ}\\ \frac{\sigma \eta \mu }{\delta ϵ}& \frac{\sigma \eta \mu }{\delta ϵ}& 0\end{array}\right)$ (9)

$V=\left(\frac{\partial \Phi \left(X\right)}{\partial X}\right)\left(E_0\right)=\left(\begin{array}{ccc}\theta +\mu & 0& 0\\ 0& \omega +\mu & 0\\ 0& 0& ϵ\end{array}\right)$

and its inverse,

$V^{-1}=\left(\begin{array}{ccc}\frac{1}{\theta +\mu }& 0& 0\\ 0& \frac{1}{\omega +\mu }& 0\\ 0& 0& \frac{1}{ϵ}\end{array}\right)$ (10)

i.e., the Next-Generation Matrix is,

$G=F\cdot V^{-1}=\left(\begin{array}{ccc}0& 0& \frac{f\beta \delta ϵ}{\mu \eta ϵ}\\ 0& 0& \frac{\left(1-f\right)\beta \delta ϵ}{\mu \eta ϵ}\\ \frac{\sigma \eta \mu }{\delta ϵ\left(\theta +\mu \right)}& \frac{\sigma \eta \mu }{\delta ϵ\left(\omega +\mu \right)}& 0\end{array}\right)$ (11)

and its corresponding characteristic equation $\left|G-\lambda I\right|=0$ , is

$\lambda \left[{\lambda }^2-\frac{f\beta \sigma }{\left(\theta +\mu \right)ϵ}-\frac{\left(1-f\right)\beta \sigma }{\left(\omega +\mu \right)ϵ}\right]=0$

which has the following eigen values, ${\lambda }_1=0$ and

${\lambda }_{1,2}=\pm \sqrt{\frac{f\beta \sigma }{\left(\theta +\mu \right)ϵ}+\frac{\left(1-f\right)\beta \sigma }{\left(\omega +\mu \right)ϵ}}$ (12)

Therefore, the spectral radius of G is,

$\rho \left(G\right)=\max \left\{0,\pm \sqrt{\frac{f\beta \sigma }{\left(\theta +\mu \right)ϵ}+\frac{\left(1-f\right)\beta \sigma }{\left(\omega +\mu \right)ϵ}}\right\}$ (13)

and the threshold, basic reproductive number,

$R_0\left(f\right)=\sqrt{\frac{f\beta \sigma }{\left(\theta +\mu \right)ϵ}+\frac{\left(1-f\right)\beta \sigma }{\left(\omega +\mu \right)ϵ}}$ (14)

i.e., $R_0\left(f\right)=\sqrt{R_0^a\left(f\right)+R_0^s\left(f\right)}$ , where

$R_0^a\left(f\right)=\frac{f\beta \sigma }{\left(\theta +\mu \right)ϵ},R_0^s\left(f\right)=\frac{\left(1-f\right)\beta \sigma }{\left(\omega +\mu \right)ϵ}$ (15)

where $R_0^a\left(f\right)$ and $R_0^s\left(f\right)$ is the number of secondary cases for asymptomatic infected people and the number of secondary cases for symptomatic infected people respectively, as function of the asymptomatic population fraction.

4. Simulation of the Model

Figure 2, shows the graphics of the lineal epidemic thresholds in $f$ for asymptomatic and symptomatic population and the basic reproduction number R₀, dependent on the asymptomatic population fraction. The system is simulated with Maple software for a 30% and 70% of asymptomatic population.

According to Figure 3, for 30% of asymptomatic population, itself reaches its maximum in 81000 individuals about 10 days later (black line), stabilizing at 90 days, while the population of symptomatic reaches a maximum value in 98000

individuals in 7 days (red line), increasing and decreasing fast to stabilize without infection in 30 days. In addition, it shows the behavior of the population of transmitting mosquitoes, for a 30% of asymptomatic population. The population grows up fast to a maximum value of 290 transmitting mosquitoes to 180 days slowly decreasing to zero in more than 100 days (blue line).

Figure 4, shows a population with a very fast growth of the asymptomatic population up to a maximum value of 195000 individuals in 10 days (black line), decreasing to stabilize in a low level after 90 days, while the symptomatic population stabilizes without infection around the 30 days (red line). For 70% of asymptomatic human population, the transmitting mosquitoes increase fast up

to reach its maximum at 350 transmitting mosquitoes around 28 days later (blue line) and from this value slowly decreases.

5. Result and Conclusions

To determine R₀, the infectious stages are considered, i.e. $x_3,x_7,x_2,x_6$ y $y_2$ . Following the approach of the Next-Generation Matrix [28] [29] [30] , we obtain the basic reproduction number, given by:

$R_0=\sqrt{\frac{\sigma {\beta }_{1}f}{ϵ\left({\alpha }_1+\mu \right)}+\frac{{\beta }_2\sigma \left(1-\left(f+g+h\right)\right)}{ϵ\left({\alpha }_2+\mu \right)}}$ (16)

To understand better the meaning of the basic reproduction number, the additive and multiplicative effects of R₀ indicate that vector can transmit the

DENV or CHIKV to the susceptible population. The terms $\frac{{\beta }_{1}f}{ϵ}$ and $\frac{{\beta }_2\left(1-\left(f+g+h\right)\right)}{ϵ}$ indicate the incidence of new cases of DENV or CHIKV in the susceptible population respectively, during the lifetime of vector. The expression $\frac{\sigma }{{\alpha }_1+\mu }$ corresponds to the DENV incidence in nontransmitting mosquitoes. On the other hand, the expression $\frac{\sigma }{{\alpha }_2+\mu }$ represents the CHIKV incidence in nontransmitting mosquitoes.

Acknowledgements

Thanks Grupo de Modelación Matemática en Epidemiología (GMME), Facultad de Educación, Universidad del Quindío―Colombia.

Cite this paper

Jamboos, J.C., Muñoz, D.M., Muñoz, A., Manrique, O.A. and Raigosa, S. (2018) Simulation Model to the Zika Virus Considering Asymptomatic Population. Open Journal of Modelling and Simulation, 6, 1-12. https://doi.org/10.4236/ojmsi.2018.61001

References