AN APPLICATION OF NONSTANDARD FINITE-DIFFERENCE SCHEME FOR SOLVING AUTONOMOUS AND NON-AUTONOMOUS MATHEMATICAL MODEL FOR WOLBACHIA-CARRYING MOSQUITO POPULATION DYNAMICS

unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract. The use of Wolbachia bacterium has been proposed as an alternative strategy against Dengue, Zika and Chikungunya. This requires that Wolbachia-carrying mosquitoes should persist in the population. A number of mathematical models has been developed and analysed to understand Wolbachia-carrying mosquito population dynamics. However, their analytical solutions are not easily derived and therefore, a numerical approach is required. In this paper, we develop a nonstandard finite difference scheme (NSFDS) for autonomous and non-autonomous mathematical models of Wolbachia-carrying mosquito population. The dynamical properties of discrete systems are then analysed. We also perform numerical simulations of the scheme and compare to other traditional methods. We found that the discrete system preserves properties of the continuous models such as equilibrium points and stability.


INTRODUCTION
The use of Wolbachia bacterium has been proposed as a promising strategy against Dengue, Chikungunya, and West-Nile virus. There are two mechanism by which Wolbachia can reduce the transmission of vector-borne diseases. First, the Wolbachia can reduce the mosquito lifespan up to 50% depending in the Wolbachia strain [1], which minimises the chance for mosquitoes to transmit virus. Second, Wolbachia can reduce the level of virus in the mosquitos [2]. Wolbachia causes an effect known as Cytoplasmic Incompatibility (CI). That is, Wolbachia-carrying female mosquitoes can reproduce successfully when mating with both non-Wolbachia and Wolbachia carrying male mosquitoes. On the other hand, non-Wolbachia female mosquitoes can only reproduce successfully when mating with non-Wolbachia male mosquitoes [2]. This gives a reproductive advantage for Wolbachia-carrying female mosquitoes. Furthermore, CI affects the population dynamics of mosquitoes.
Mathematical model is commonly formulated and analysed to understand complex phenomena including population dynamics [3,4,5,6,7,8,9,10]. The models are generally in the form of system of nonlinear differential equations, which consist of autonomous or non-autonomous models. Furthermore, analytical solutions of the models are not easily determined and therefore, a numerical approach is generally used [11]. To date, a number of mathematical models has been developed to understand the mosquito population dynamics in the presence of Wolbachia-carrying mosquitoes [4,12,3]. The model are autonomous and non-autonomous with seasonal forcing on adult mosquito death rate [4]. Analytical solutions of the models are not easily determined and therefore numerical solutions are presented. Although, NSFDS has been formulated for other problems [13,14,11,15,16,17,18,19,20,21,22,23], to the best of our knowledge, there is a little work on the construction of NSFDS for model of Wolbachia-carrying mosquito population dynamics. The aim of this paper is to construct a NSFDS for autonomous and nonautonomous mathematical model of Wolbachia-carrying mosquito population and analyse the scheme properties.
The remainder of the paper is organised as follows. Section 2 overviews mathematical model for Wolbachia-carrying mosquito population dynamics. Section 3 presents the construction of numerical scheme for autonomous and nonautonomous models and their numerical simulations.
Finally, the conclusion is presented.

MATHEMATICAL MODEL
In this section, we present a mathematical model for Wolbachia-carrying mosquito population dynamics. The model has been proposed by Ndii et al. [4]. The mosquito population consists of non Wolbachia and Wolbachia-carrying mosquitoes. The mosquito population is divided into Aquatic (A N and A W ), Male (M N and M W ) and female (F N and F W ) mosquitoes.
The subscripts N and W are to differentiate between non-Wolbachia and Wolbachia-carrying mosquitoes respectively.
The Wolbachia-carrying female mosquitoes reproduce when they mate with non-Wolbachia and Wolbachia-carrying male mosquitoes and their growth is limited by carrying capacity, K, that is Non-Wolbachia female mosquitoes can only reproduce when they mate with non-Wolbachia males. This is governed by Equations (1) and (2) capture the effect of Cytoplasmic Incompatibility (CI). Furthermore, the maternal transmission of Wolbachia-carrying mosquitoes is not perfect [2] and hence there is a proportion of Wolbachia-carrying aquatic mosquitoes (1 − α) that mature to be non-Wolbachia mosquitoes. A proportion of α mature to be Wolbachia-carrying adult mosquitoes. When aquatic mosquitoes mature, a proportion of ε is male and the rest (1 − ε) is female. The model is governed by the following system of differential equations.
where P = F N + M N + M W + F W . We nondimensionalise the model by setting the ratio of male and female mosquitoes to be the same, ε = 1/2 and K = 1, the model is reduced to In this paper, a NSFDS is constructed for nondimensionalised model (Equation (4)).
Model (4) has four equilibriums, which are mosquito-free equilibrium, only non-Wolbachia mosquito surviving equilibrium, only and non-Wolbachia and two coexistence equilibriums.
The expression for E 3,4 is not analytically tractable and hence is explored numerically as given in Ndii et al. [4]. Furthermore, only one of them is locally stable.
As the mosquito lifespan is seasonal-dependent, a seasonal forcing is applied to the adult mosquito death rate. Furthermore, a sensitivity analysis showed that adult mosquito death rate is also the influential parameter [4]. The seasonal forcing of the mosquito adult death rate is the following µ N (t) = µ N0 (1 + η cos(2π(t + ψ)) and µ W (t) = rµ N (t) µ N0 is the baseline adult death rate, η is the degree of seasonality, ψ is the phase and t is time, and r is the ratio of death rate of Wolbachia-carrying adult mosquitoes to non-Wolbachia adult mosquitoes.

SCHEME CONSTRUCTION
This section deals with numerical construction of the proposed scheme. A nonstandard numerical scheme is based on the two fundamental rules [27,28,29], which are (1) Nonlocal approximation is used; for example (2) Discretisation of derivatives is not traditional and use the nonnegative function φ (h) = Let us define the derivatives as follows is real-valued function on ℜ called denominator function that satisfies the following properties [30] ( The above definition is consistent with the traditional derivatives. There is no general rule for determining φ (h) but we can find ideas in [30].
Rearrange the Equation (6) to obtain We choose the denominator function as

Scheme analysis.
This section presents the equilibrium points of the scheme and the stability of the model. We determine the equilibrium points that satisfy the conditions of A n+1 , If Wolbachia-carrying mosquitoes do not persist in the population, then the equilibrium points is The equilibrium E 2 corresponds exactly to that of Model (4).
The other equilibrium points (E 3 and E 4 ) can be found numerically since the mathematical expressions are not analytically tractable.
To determine the stability, we construct a Jacobian matrix and find the eigenvalues. For discrete model, the equilibrium points are stable if the root of characteristic polynomial is less than unity. Constructing the Jacobian matrix J(x), which is , The eigenvalues of the model is found numerically by solving |J − λ I|. The results are given in Table 2.  Table 2 . We show the results that the absolute of the eigenvalues is always less than unity and therefore the equilibrium points are locally stable. Note that the obtained equations should be computed in sequence but it starts from F n+1 for the calculation of A n+1 N , which is then used to calculate A n+1 W and then F n+1 W . The numerical simulations are given in Section 3.3

Numerical simulations.
To illustrate the results, we perform numerical simulations of the scheme and compare to the Euler method. The parameter values are given in Table 1     where r is the ratio of death rate of Wolbachia-carrying adult mosquitoes to that of non-Wolbachia adult mosquitoes. The numerical scheme for the nonautonomous model is the following.
where their expressions are the same as Equation (7) except for µ N and µ W . The µ N and µ W are replaced by Equations 10. This gives the periodic solutions of the model. The simulations of the model is given in Section 3.5. The calculation of the obtained equations begins by calculating the parameter µ n+1 N . The following steps are the same as in autonomous model.

Numerical Simulations.
In this section, we present numerical simulations of NSFDS for the nonautonomous model and compare the results to that of the Euler method and matlab ode45 routine. The simulation is given in Figure 3.

CONCLUSIONS
We formulated a nonstandard finite difference schemes for autonomous and non-autonomous model of Wolbachia-carrying mosquito population. A seasonal forcing is added in adult mosquito death rate as this is influential parameter of the model. We found that the scheme is locally stable and the numerical simulations are in good agreement with Euler method and Matlab ode45 routine. For large time window, the numerical simulations show that the solutions remain positive as long as φ i (h) < 1. However, the convergence of the scheme depends on the step size h. Furthermore, one may be interested to use other denominator functions such as hyperbolic tangent function. This may improve the accuracy of the scheme.