Reaction-diffusion Model Applied to the Population Dynamics of Wild and Transgenic Mosquitoes

Due to recent advances in genetic manipulation, transgenic mosquitoes can be a viable alternative to reduce some diseases. Viability conditions are obtained by the simulation and analysis of mathematical models that describe the behavior of wild and transgenic mosquitoes population living in the same geographic area. In this work, we present a reaction-diffusion model with a nonlinear reaction term, a function that describes the interaction between wild and transgenic mosquitoes taking into account their zygosity. The diffusive term represents a uniform spatial spread characterized by a fixed diffusion parameter. The system of partial differential equations obtained is solved numerically by combining a implicit Runge-Kutta method and finite elements method, through the sequential operator splitting technique. Several scenarios are analyzed, simulating the spatial release of transgenic mosquitoes, and lead to an understanding of an intrinsic relationship between the transgenic and wild varieties for different initial conditions.


INTRODUCTION
Vector-borne diseases have always been a big preoccupation for populations and government authorities in tropical countries, especially in those with low human development rates. The combination of a favorable climate and intensive agriculture, associated with deforestation and poor sanitation conditions, provide ideal environmental conditions for the proliferation of many species of vectors responsible for the transmission of several diseases that affect the people of Genetically modified mosquitoes refractory to malaria were first obtained in 2002, using a technique developed by Catteruccia et al. [1]. To obtain them, the scientists developed two different types of Anopheles stephensi using the CP (carboxypeptidase) promoter: one of them expressing peptide SM1 (salivary gland and midgut binding peptide 1) [10] and the other expressing the enzyme PLA2 (phospholipase A2), present in bee venom [16]. These new insects must interact with wild mosquitoes by mating and spreading the gene that determines the interruption of the transmission process.
In 2015, Gants et al. and Bier et al. [5] developed a method based on the self-propagating CRISPR/Cas9 genome-editing technology that converts heterozygous mutations to homozygotes. This mechanism of artificial activation, called mutagenic chain reaction (MCR), was tested for Drosophila melanogaster with a satisfactory level of 97%. Gantz et al. [6] and Noble et al. [18] proved to be efficient for A. stephensi, although results obtained by mating between transgenic males and wild-type females were more efficient than mating results between wild males and transgenic females.
The modeling of the dynamics conducting the interaction between wild and transgenic mosquitoes has been studied ever since. Mathematical and computational models can be found in the literature, highlighting the works developed by Li et al. [12], Diaz et al. [8] and Wyse et al. [24]. In 2004, Li et al. [12] presenting a discrete model that considered the interaction between two varieties of mosquitoes: wild and transgenic, without distinction of zygosity; in this model the genetic mutation was considered without reduction or favoritism in the vital rates of the mosquitoes. The rates of birth and mortality were considered dependent of density and two situations were analyzed: constant mating rate and mating rate proportional to the total population. A continuous-time version, described by a system of ordinary differential equations, was obtained by Li et al. in [13]. A functional response type Holling II and the Alee effect in wild The fitness of wild and transgenic mosquitoes was considered in the model of Diaz et al. [8]. This model, in continuous time, also took into account the zygosity of transgenic mosquitoes. A discrete time model, without distinction of zygosity, taking into account the horizontal and vertical transmission of a genetically modified bacterium was proposed by Li in [14], assuming horizontal transmission depends on the mating between wild mosquitoes and those that have mutated due to contact with the bacteria.
In this paper, a mathematical model that describes the interaction dynamics through mating between transgenic and wild mosquitoes is investigated, as well as the dissemination of the transgene that determines the interruption of an epidemiological process. For this, the transgenic mosquitoes are differentiated according to their zygosity, being heterozygous or homozygous. The interaction between mosquitoes describes a density dependency for vital rates and imposes a maximum limit of population growth that occurs from the carrying capacity. The diffusion is represented by Fick's law, with the fixed diffusion coefficient. Thus, the mathematical model obtained is a non-linear diffusion-reaction system. The resulting system is solved numerically by the operator splitting method, which is well known in solving problems resulting in large systems of partial differential equations, as well as problems involving nonlinear chemical reactions [2], mosquito dispersion [3,24] and non-linear applications of Schrodinger [19].
In section 2, the mathematical model that describes the population dynamics resulting from the introduction of genetically modified mosquitoes in wild populations is presented. The proposed mathematical model is based on strategies aimed at genetically modified mosquitoes that are designed to have a reduced transmission capacity of a given infectious agent. They are also fertile and able to propagate and perpetuate their hereditary trait in the wild mosquito population. In section 3, the methodology used to obtain the numerical solution is shown, and consists in the technique of decomposition of operators. The fourth order Runge-Kutta method is used for the dynamic problem and the finite element method is used for the spatial problem. Section 4 contemplates the numerical results obtained from different initial conditions. Finally, in section 5 there is a brief discussion and some suggestions for future investigations are pointed out.

MATHEMATICAL MODELING OF MOSQUITO DISPERSAL
The mathematical model presented is based on the Fisher-KPP equation [4,11], assuming that the populations have the same fitness, according to studies on A. stephensi [17], which ensures to the mosquitoes the reproductive success and adaptations to the environment wild, with competitive ability. Each individual progresses and emerges in adulthood at a net rate and leaves due to mortality. To obtain the reaction terms, we consider the population dynamics of mosquitoes governed by the logistic equation with capture: where N is the total population of adult females, γ = ε − δ 1 is the net inflow rate of mosquitoes into adulthood, where ε is the entry rate and δ 1 is the death rate. In addition to δ 1 , there is also a mortality rate δ 2 , independent of density, introduced in order to take into account other factors inducing mortality that keeps the population stabilized at a level below the support capacity C . Thus, equation (2.1) can be written in the equivalent form: Assuming that the total population is composed of wild mosquitoes (u 1 ), heterozygous transgenics (u 2 ) and homozygous transgenics (u 3 ), so that N = ∑ u i , and holds: As mosquitoes are diploid organisms, it can establish the genotype of wild mosquitoes as (w, w) and the genotype of the homozygous transgenic mosquitoes as (g, g). Denote by a i j , b i j and c i j the genotypic frequencies for u 1 , u 2 and u 3 obtained from mating u i × u j , i, j = 1, 2, 3. These coefficients satisfy the relation a i j + b i j + c i j = 1.
Given a i j , b i j , c i j , ε, γ,C, δ , κ ∈ R, the problem is to find u i (x,t) ∈ R for all x ∈ Ω and t > t 0 ∈ R + through the resolution of the following reaction-diffusion system: with Dirichlet contour conditions given by where t 0 ∈ R + is the initial time instant and κ is the fixed diffusion coefficient since there is no evidence of flight alteration caused by genetic manipulation of this species.

NUMERICAL FORMULATION
In this section, the focus is on the development of discrete formulations and applications of computational techniques to numerically solve the proposed model. For this, it is used a technique of sequential operator splitting [7] to dissociate the original system into another equivalent, formed by a combination of two subsystems that fall into problems of less complexity. Taking advantage of the fact that the structure of the reaction-diffusion models allows a natural decomposition of the equations, thus providing the opportunity to apply schemes of operator splitting. Thus, it becomes possible to separately handle each of the systems and solve each problem with the numerical method that best fits the nature of the operator involved.
To describe the algorithm, it is necessary to proceed with the decomposition of the system (2.4) into two problems: a system of partial differential equations with the purely diffusive term, and the other of nonlinear ordinary differential equations, associated to the purely reactive term. To solve the diffusive problem let us use the finite element method with an implicit finite difference scheme of the type Crank-Nicolson and Euler. To solve the second system, which is nonlinear, let us the fourth-order Runge-Kutta method.
Let us denote Ω ≡ (0, L) the spatial domain and I ≡ (0, T ) the time interval of interest, with T > 0 the final time. By introducing the time discretization [0, T ] = N n=0 [t n ,t n+1 ], with I t n = [t n ,t n+1 ] a partition of I,N = T /∆t the number of partitions and ∆t = t n+1 −t n a time step, which we consider uniform. With this, follow the steps: Step 1: For the initial time, t = t 0 , initialize the variables u i (x,t 0 ) = g i (x), for each i = 1, 2, 3, being g i (x) the initial condition given.
Step 2: For a fixed n, n ≥ 0, given the initial conditions u i (x,t n ), calculate u i (x,t) at time t n+1 through the following problem: Problem A: Given D ∈ R, find u i (x,t) ∈ R, with x ∈ Ω and t ∈ I n satisfying the system: with boundary conditions and initial conditions Step 3: In the same time interval I t n , use the solution of Problem A, given by the previous step, as the initial condition for calculating the solution of the system of coupled nonlinear ordinary differential equations associated with the model (2.4), expressed by the following problem: with initial conditions where u i (x,t n+1 ) are the solutions obtained from Problem A.
Step 4: The solution of Problem B is the approximate solution of the model at time t n+1 ∈ I t n ⊂ I. If t n+1 < T , increments n, returns to Step 2 and repeats the process until equality occurs.
For the resolution of Problem A, a finite element method approach is used. For this, consider a Sobolev space and a Hilbert space, respectively. The spatial domain is discretized using a uniform finite element partition consisting of n e elements Ω e , such that Ω = n e n=1 Ω e and n e n=1 Ω e = / 0. This choice is to construct are Lagrange polynomials.
Consider the Problem (A h ): Let u h i (x,t) ∈ U h , (i = 1, 2, 3), t ∈ I n , such that: with initial conditions given by To solve this problem, a transient algorithm is required to obtain the numerical solution of the semi-discrete problem. For this, take the approximation at a time t as follows: and v h (x,t) = ∑ n p j=1 ϕ j (x) to construct the interpolation functions in finite elements, respectively, where n p is the total number of freedom degrees and ϕ j (x) are the global form functions. Thus, the problem (A h ) leads to the following system of ordinary differential equations: where which can be written in the following matrix form as: To solve this system numerically, it is enough to use the transient algorithm, based on the generalized family of trapezoidal methods [9]: for α = 0, α = 0, 5 or α = 1 it has, respectively, the methods of Forward Euler, Crank-Nicolson or Backward Euler.
In relation to Problem B given by the system of ordinary differential equations (3.4), use the fourth order Runge-Kutta method.

NUMERICAL SIMULATIONS
In this section, the obtained results are reported when applying computational techniques to numerically solve the proposed reaction-diffusion model. It is an analysis of the behavior of the three mosquito varieties in different scenarios, considering Mendelian (as SM1 mutation [16])).
Three situations are considered for numerical experiments, as well as population dynamics with different initial conditions: • The first experiment attempts to establish a correlation with a controlled laboratory experiment described by Moreira et al. [16]. It considers an identical initial quantitative for the population of wild mosquitoes and heterozygous transgenic mosquitoes, releasing heterozygous transgenic mosquitoes closer to the wild mosquitoes.
• The second experiment establish different mosquitoes release scenarios to establish the most effective release strategies. It considers the release of mosquitoes into two distinct regions of the domain. In practice, this idea consists of identifying the main outbreaks of wild mosquitoes, locating them on the integration domain and inserting transgenic mosquitoes into these identified regions.
• The third experiment considers an initial conditions analogous to the first experiment, but with the initial population composed only of wild mosquitoes and homozygous transgenic. The population of heterozygous mosquitoes will initially be null, and it will emerge from the mating process among the varieties considered.
In all experiments, it was considered that mosquitoes occupy a spatial region I d = [0, 10], measured in km, whose propagation starts from the release of mosquitoes and it has a duration of four weeks, I t = [0, 4]. These intervals are large enough for populations to diffuse without reaching the extremes of this range, which would lead to a spontaneous loss of mosquitoes as a function of zero contour conditions.
All the parameters used in the numerical simulations are in Tables 2 and 3.   It was arbitrarily adopted the diffusion coefficient κ, based on the fact that the flight range of the mosquitoes is small, since Anopheles spp. acts in the intra and peridomestic. As there are no reports that indicate the effects of genetic manipulation on the flight mode of mosquitoes, so it was considered the same value of κ for the three varieties considered.

Release of heterozygous into a wild mosquito focus
In these simulations will be considered an initial quantitative identical for the population of wild mosquitoes and heterozygous transgenic mosquitoes, distributed over the range Ω = [0, 10] in two distinct ways: in the first, it will simulate the behavior of populations when it is released heterozygous transgenic mosquitoes in a focus of wild mosquitoes; in the second it will simulate this behavior when it is released heterozygous transgenic mosquitoes in a region where the wild population is homogeneously distributed. Initially, it will verify the dynamic behaviour of the three interacting populations. In Figure 1, the initial populations of equally distributed heterozygous wild and transgenic mosquitoes were considered and the initial population of homozygous transgenes was equal to zero. These conditions correspond to the laboratory test performed by Moreira et al. [16], resulting in a stabilization of the population with 56% of wild mosquitoes and 44% of transgenic mosquitoes, ratified by the Hardy -Weinberg equilibrium. In this simulation, no transgenic homozygous individual was released, but it arises naturally from the mating process. The initial conditions used for the wild, heterozygous and homozygous transgenic populations were u 1 (0) = u 2 (0) ≈ 0, 8 and u 3 (0) = 0, obtained from (4.1) and (4.2). In Figure 2, the system (2.4) was simulated satisfying Dirichlet boundary conditions and initial conditions (4.3) and (4.4). These conditions are quantitatively equivalent to those used in the simulation of Figure 1. However, this quantity was distributed along the spatial domain. This situation represents the introduction of heterozygous transgenic mosquitoes into a wild mosquito focus, located at the center of the integration domain.

Release of wild and transgenic heterozygous mosquitoes in two positions
In this section, the release of mosquitoes into two distinct regions of the Ω domain it will be considered. In practice, this idea consists of identifying the main outbreaks of wild mosquitoes, locating them on the integration domain and inserting transgenic mosquitoes into these identified regions.
In Figure 3, initial populations of wild mosquitoes and heterozygous transgenic mosquitoes of the same size and spatial distribution, both with higher density in the positions x = 4.75 and x = 6.25 were considered. In field conditions, this could be done by identifying two regions with the highest incidence of wild mosquitoes, locating this region in the integration interval and releasing an equivalent amount of heterozygous transgenic mosquitoes at this local.

Release of wild mosquitoes and homozygous transgenics
In this section, it will consider the initial population composed only of wild mosquitoes and homozygous transgenic. The population of heterozygous mosquitoes will initially be null, but it will emerge from the mating process among the varieties considered. All simulations presented in this section will have initial conditions analogous to those in subsection 4.1, only inserting homozygous transgenic instead of heterozygous. In this simulation the curves that representing the populations of wild and transgenic homozygous mosquitoes are superimposed. It is not actually possible to see it that populations of wild and transgenic homozygotes mosquitoes stabilized at the same level. Corresponding to 25% of the total population of wild mosquitoes, 25% of homozygous transgenic mosquitoes and 50% of heterozygous transgenic mosquitoes.
In Figure 5 the system (2.4) was simulated satisfying Dirichlet boundary conditions and initial conditions 4.9) and (4.10). These conditions represent an introduction of homozygous transgenic mosquitoes into a focus of wild mosquitoes, with the highest concentration located at the center of the integration range. The amount of homozygous transgenic mosquitoes introduced is equal to that observed for wild mosquitoes.
The value obtained for the density of each population at the end of integration, obtained by the trapezoid rule, is u 1 (4) = 19.6735, u 2 (4) = 39.3471 and u 3 (4) = 19.6735. The proportions obtained are consistent with those obtained in Figure 4.

CONCLUSION
In this article, the development of a mathematical model that describes the interaction between wild and transgenic mosquitoes is presented, taking into account zygosity, considering all populations in absolute numbers and admitting a fixed diffusion coefficient. These characteristics make the description of the model more realistic, ensuring the existence of all possible phenotypes, preserving the parameters that do not exist in the dimensionless form and allowing to the mosquitoes the same vital rates and capacity of uniform displacement.
The dynamic system was developed to preserve the peculiarities of the species and to avoid overlapping of individuals when the transgenics are inserted, respecting the support capacity of the environment. The terms of the dynamic system related to the mating and competition for resources implied in the nonlinearity of this system, which is characteristic of the great majority of population dynamics models. The diffusion term, based on Fick's law, describes a symmetrical spread of the mosquito's population as a Gaussian process.
The solution of the proposed problem was obtained using the operator splitting method to decouple the reaction-diffusion system into two subproblems. The diffusive problem was solved using the finite element method and the reactive problem using the Runge -Kutta method fourth order. The low computational cost, easy of implementation and the history of good results already obtained in the modeling of chemical reactions encouraged us to adopt this procedure.
The results shown are consistent with the expected behavior, indicating a constant presence of transgenic mosquitoes after entering into the ecosystem and reducing the costs of periodic insertion. Without the superiority of transgenics, the total elimination of wild mosquitoes is impossible to achieve, since they are also obtained from the mating between the heterozygous transgenics. It is worth mentioning that the extinction of a species should not be the main objective, it is sufficient that the population of wild mosquitoes is reduced to levels that are not harmful to human health.
The proposed model with adequate epidemiological dynamics is able to study the impact of the reduction of infectious diseases transmitted by mosquitoes. By identifying an effective breeding site, it is possible to devise an optimal strategy for the release of transgenic mosquitoes.
This is an initial study that opens possibilities for future investigations. Initial conditions more adequate to reality need to be tested to obtain a strategy for the release of transgenic insects in the environment, by means of genetic algorithms, for example, in mono and multiobjective versions, considering factors such as seasonality, temperature, and rainfall.