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

ABSTRACT

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.

Keywords:
mathematical model; operator splitting; genetically modified mosquito

RESUMO

Devido aos progressos recentes na manipulação genética, os mosquitos transgênicos podem ser uma alternativa viável para reduzir algumas doenças. Condições de viabilidade são obtidas pela simulação e análise de modelos matemáticos que descrevem o comportamento de populações de mosquitos selvagens e transgênicos convivendo em uma mesma área geográfica. Neste trabalho, apresentamos um modelo de reação-difusão com um termo de reação não linear, uma função que descreve a interação entre mosquitos selvagens e transgênicos, levando em consideração sua zigosidade. O termo difusivo representa uma dispersão espacial uniforme caracterizada por um parâmetro de difusão fixo. O sistema de equações diferenciais parciais obtidas é resolvido numericamente combinando um método implícito de Runge-Kutta e um método de elementos finitos, através da técnica de divisão sequencial de operadores. Vários cenários são analisados, simulando a liberação espacial de mosquitos transgênicos, e levam ao entendimento de uma relação intrínseca entre as variedades transgênica e selvagem para diferentes condições iniciais.

Palavras-chave:
modelo matemático; decomposição de operadores; mosquitos geneticamente modificados

1 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 these regions. Due to the acceleration of global warming, especially in the last decade, this concern has extended to regions beyond the tropics, making many of these diseases a focus of world attention ²⁰.

According to the World Health Organization ²², about 17% of the global population suffers from infectious diseases transmitted by mosquitoes, malaria being the most known and deadly disease. In 2016, 91 countries reported a total of 216 million malaria cases, resulting in 445.000 deaths, approximately the same number reported in 2015 with 395.000 deaths ²¹. However, among vector-borne diseases, the incidence of dengue is the fastest growing in the world, with a 30-fold increase in the last 50 years. In addition to malaria and dengue, yellow fever, chikungunya, zika and others are also stand out.

Recent advances in the genetic manipulation of mosquitoes, in particular, Aedes spp. and Anopheles ssp., aim at new approaches to control vector-borne diseases. In this way, the use of genetically modified mosquitoes can be a viable alternative for the control of mosquito-borne infections, acting together with prophylaxis, vaccines, insecticides and medicines.

Genetically modified mosquitoes refractory to malaria were first obtained in 2002, using a technique developed by Catteruccia et al. ¹. 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) ¹⁰ and the other expressing the enzyme PLA2 (phospholipase A2), present in bee venom ¹⁶. 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. ⁵ 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. ⁶ and Noble et al. ¹⁸ 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. ¹², Diaz et al. ⁸ and Wyse et al. ²⁴. In 2004, Li et al. ¹² 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 ¹³. A functional response type Holling II and the Alee effect in wild and transgenic populations was considered in this model. A continuous time model with zygosity distinction was also proposed by Wyse et al. in ²³.

The fitness of wild and transgenic mosquitoes was considered in the model of Diaz et al. ⁸. 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 ¹⁴, 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 ², mosquito dispersion ^{3), (24} and non-linear applications of Schrodinger ¹⁹.

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.

2 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¹⁷, 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:

$\frac{dN}{dt}=\gamma N\left(1-\frac{N}{C}\right)-{\delta }_{2}N,$ (2.1)

where N is the total population of adult females, $\gamma =\epsilon -{\delta }_1$ is the net inflow rate of mosquitoes into adulthood, where ε is the entry rate and δ1 is the death rate. In addition to δ ₁, there is also a mortality rate δ ₂, 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:

$\frac{dN}{dt}=\epsilon N-\left(\frac{\gamma }{C}\right)N^2-({\delta }_1+{\delta }_2)N=\left(\frac{\epsilon }{N}-\frac{\gamma }{C}\right)N^2-\delta N,$ (2.2)

where $\delta ={\delta }_1+{\delta }_2$.

Assuming that the total population is composed of wild mosquitoes (u ₁), heterozygous transgenics (u ₂) and homozygous transgenics (u ₃), so that $N=\sum u_i$, and holds:

$\frac{d(\sum u_i)}{dt}=\left(\frac{\epsilon }{\sum u_i}-\frac{\gamma }{C}\right){\left(\sum u_i\right)}^2-\delta \sum u_i.$ (2.3)

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 _ij , b _ij and c _ij the genotypic frequencies for u ₁, u ₂ and u ₃ obtained from mating $u_i\times u_j,i,j=1,2,3$. These coefficients satisfy the relation $a_{ij}+b_{ij}+c_{ij}=1$.

The Table 1 synthesizes the wild-type mosquito populations (w, w), homozygous transgenic (g, g) and heterozygous (w, g) after the crosses.

Table 1:
Genotypic frequencies obtained from crosses between wild, heterozygous and homozygous mosquitoes.

Crossing	$u_1\times u_1$	$u_1\times u_2$	$u_1\times u_3$	$u_2\times u_2$	$u_2\times u_3$	$u_3\times u_3$
Genotype
(w,w)	a 11	a 12	a 13	a 22	a 23	a 33
(w,g)	b 11	b 12	b 13	b 22	b 23	b 33
(g,g)	c 11	c 12	c 13	c 22	c 23	c 33

Given $a_{ij},b_{ij},c_{ij},\epsilon ,\gamma ,C,\delta ,\kappa \in \mathrm{ℝ}$, the problem is to find $u_i(x,t)\in \mathrm{ℝ}\text{for all}x\in \Omega$ and $t>t_0\in {\mathrm{ℝ}}^{+}$ through the resolution of the following reaction-diffusion system:

$$\begin{gathered} \left\{\frac{\partial u_1}{\partial t}=\kappa \frac{{\partial }^2{u}_1}{\partial x^2}+\left(\frac{ϵ}{\sum u_i}-\frac{\gamma }{C}\right)\sum \sum a_{ij}{u}_i{u}_j-\delta u_1, \\ \frac{\partial u_2}{\partial t}=\kappa \frac{{\partial }^2{u}_2}{\partial x^2}+\left(\frac{ϵ}{\sum u_i}-\frac{\gamma }{C}\right)\sum \sum b_{ij}{u}_i{u}_j-\delta u_2, \\ \frac{\partial u_3}{\partial t}=\kappa \frac{{\partial }^2{u}_3}{\partial x^2}+\left(\frac{ϵ}{\sum u_i}-\frac{\gamma }{C}\right)\sum \sum c_{ij}{u}_i{u}_j-\delta u_3,\right. \end{gathered}$$ (2.4)

with Dirichlet contour conditions given by

$u_i(x_0,t)={\overline{u}}_{i0},u_i(x_L,t)={\overline{u}}_{iL},$ (2.5)

and initial conditions

$u_i(x,t_0)={\overline{u}}_i\left(x\right),$ (2.6)

where $t_0\in {\mathrm{ℝ}}^{+}$ 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.

3 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 ⁷ 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 $\Omega \equiv (0,L)$ the spatial domain and $I\equiv (0,T)$ the time interval of interest, with $T>0$ the final time. By introducing the time discretization $[0,T]={\cup }_{n=0}^N[t_n,t_{n+1}]$, with $I_{t_n}=[t_n,t_{n+1}]$ a partition of $I,N=T/\Delta t$ the number of partitions and $\Delta 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 ${\overline{u}}_i(x,t_0)=g_i\left(x\right)$, for each $i=1,2,3$, being $g_i\left(x\right)$ the initial condition given.

Step 2: For a fixed n, $n\ge 0$, given the initial conditions ${\tilde{u}}_i(x,t_n)$, calculate ${\tilde{u}}_i(x,t)$ at time $t_{n+1}$ through the following problem:

Problem A: Given $D\in \mathrm{ℝ}$, find ${\tilde{u}}_i(x,t)\in \mathrm{ℝ}$, with $x\in \Omega$ and $t\in I_n$ satisfying the system:

$\frac{\partial \widetilde{u_i}(x,t)}{\partial t}=\kappa \frac{{\partial }^2\widetilde{u_i}(x,t)}{{\partial x}^2},$ (3.1)

with boundary conditions

$\widetilde{u_i}(0,t_n)={\overline{u}}_{i0},\widetilde{u_i}(L,t_n)={\overline{u}}_{iL},$ (3.2)

and initial conditions

$\widetilde{u_i}(x,t_n)=\overline{u_i}\left(t_n\right).$ (3.3)

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:

Problem B: Given $\{a_{ij},b_{ij},c_{ij},\gamma ,C,\epsilon ,\delta \}\in \mathrm{ℝ}$, find $u_i\left(t\right)\in \mathrm{ℝ}$, with $t\in I_{t_n}$, satisfying the system:

$$\begin{gathered} \left\{\frac{d{u}_1}{dt}=\left(\frac{ϵ}{\sum u_i}-\frac{\gamma }{C}\right)\sum \sum a_{ij}{u}_i{u}_j-\delta u_1, \\ \frac{d{u}_2}{dt}=\left(\frac{ϵ}{\sum u_i}-\frac{\gamma }{C}\right)\sum \sum b_{ij}{u}_i{u}_j-\delta u_2, \\ \frac{d{u}_3}{dt}=\left(\frac{ϵ}{\sum u_i}-\frac{\gamma }{C}\right)\sum \sum c_{ij}{u}_i{u}_j-\delta u_3,\right. \end{gathered}$$ (3.4)

with initial conditions

$u_i\left(t_n\right)={\tilde{u}}_i(x,t_{n+1}),$ (3.5)

where $\widetilde{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}\in I_{t_n}\subset 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 $U=\left\{u_i\right(x,t)\in H^1(\Omega )\mid u_i(0,t)={\overline{u}}_{i0};u_i(L,t)­={\overline{u}}_{iL}\}$ and $V=H_0^1\left(\Omega \right)$, 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 $\overline{\Omega }={\cup }_{n=1}^{n_e}{\Omega }_e$ and ${\cup }_{n=1}^{n_e}{\Omega }_e=\varnothing$. This choice is to construct $U_h=\left\{u_h\right(x,t)\in C^0(\Omega )\mid u_h(0)=u_h(L)=0,\forall u_h({\Omega }_e)\in \hat{P}({\Omega }_e\left)\right\}\subset U$ and $V_h\subset V$, where $\hat{P}\left({\Omega }_e\right)$ are Lagrange polynomials.

Consider the Problem (A _h ): Let ${\tilde{u}}_i^h(x,t)\in U_h,(i=1,2,3),t\in I_n$, such that:

${\int }_0^L\frac{\partial {\tilde{u}}_i^h}{\partial t}{v}_{h}dx+D{\int }_0^L\frac{\partial {\tilde{u}}_i^h}{\partial x}\frac{\partial v_h}{\partial x}dx=0,\forall v_h\in V_h,$ (3.6)

with initial conditions given by

${\tilde{u}}_i^h(x,t_n)={\hat{u}}_i\left(t_{n+1}\right).$ (3.7)

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: ${\tilde{u}}_i^h(x,t)=\sum _{j=1}^{n_p}\left({\tilde{u}}_i{)}_j{\phi }_j\right(x)$ and $v_h(x,t)=\sum _{j=1}^{n_p}{\phi }_j\left(x\right)$ to construct the interpolation functions in finite elements, respectively, where n _p is the total number of freedom degrees and ${\phi }_j\left(x\right)$ are the global form functions. Thus, the problem (A _h ) leads to the following system of ordinary differential equations:

$$\begin{gathered} \sum _{j=1}^{n_p}\left[M_{ij}\frac{d{\tilde{u}}_j}{dt}+K_{ij}\tilde{u}\right]=0,t\in I;1\le i\le n_p, \\ {\tilde{u}}_i\left(0\right)={\hat{u}}_i\left(x\right), \end{gathered}$$ (3.8)

where

$M_{ij}={\int }_{{\Omega }_e}{\phi }_i\left(x\right){\phi }_j\left(x\right)d{\Omega }_e,$ (3.9)

$K_{ij}={\int }_{{\Omega }_e}D\phi {\text{'}}_i\left(x\right)\phi {\text{'}}_j\left(x\right)d{\Omega }_e,$ (3.10)

which can be written in the following matrix form as:

$$\begin{gathered} \mathbf{\text{M}}\dot{\widetilde{\mathbf{\text{u}}}}\left(t\right)+\mathbf{\text{K}}\widetilde{\mathbf{\text{u}}}\left(t\right)=\mathbf{\text{0}}, \\ \widetilde{\mathbf{\text{u}}}\left(t_n\right)=\widehat{\mathbf{\text{u}}}\left(t_{n+1}\right). \end{gathered}$$ (3.11)

To solve this system numerically, it is enough to use the transient algorithm, based on the generalized family of trapezoidal methods ⁹:

$(\mathbf{\text{M}}+\alpha \Delta t\mathbf{\text{K}}){\mathbf{\text{u}}}_{n+1}=(\mathbf{\text{M}}-(1-\alpha \left)\Delta t\mathbf{\text{K}}\right){\mathbf{\text{u}}}_n,$ (3.12)

for $\alpha =0,\alpha =0,5$ or $\alpha =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.

4 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 ¹⁶)).

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. ¹⁶. 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 ${\overline{I}}_d=[0,10]$, measured in km, whose propagation starts from the release of mosquitoes and it has a duration of four weeks, ${\overline{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.

Table 2:
Model parameters estimated from literature data ^{15), (25}.

Parameter	Description	Value
κ	Coefficient of diffusion (km2/week)	5.10-3
ε	Emergence rate to the adult stage	8.11
δ 1	Density-dependent mortality rate	4.0
δ 2	Density-independent mortality rate	2.54
C	Support capacity	262.0

Table 3:
Genotypic frequencies considering Mendelian inheritance.

Frequency	a 11	a 12	a 13	a 22	a 23	a 33
Mendelian	1.0	0.5	0.0	0.25	0.0	0.0
Frequency	b 11	b 12	b 13	b 22	b 23	b 33
Mendelian	0.0	0.5	1.0	0.5	0.5	0.0
Frequency	c 11	c 12	c 13	c 22	c 23	c 33
Mendelian	0.0	0.0	0.0	0.25	0.5	1.0

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.

4.1 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 $\Omega =[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.

$u_1\left(0\right)=u_2\left(0\right)=\left\{\begin{array}{r}{\int }_{4.5}^{5.5}10{\left(\sin \left(\pi \left(\frac{x-4.5}{5.5-4.5}\right)\right)\right)}^{100}\text{if}4.5\le x\le 5.5,\\ 0\text{if}0\le x\le 4.5\text{and}5.5\le x\le 10.\end{array}\right.$ (4.1)

$u_3\left(0\right)=0.$ (4.2)

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. ¹⁶, 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\left(0\right)=u_2\left(0\right)\approx 0,8$ and $u_3\left(0\right)=0$, obtained from (4.1) and (4.2).

Figure 1:
Population of wild mosquitoes, heterozygous and homozygous transgenic with initial conditions (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.

$u_1(x,0)=u_2(x,0)=\left\{\begin{array}{r}10{\left(\sin \left(\pi \left(\frac{x-4.5}{5.5-4.5}\right)\right)\right)}^{100}\text{if}4.5\le x\le 5.5,\\ 0\text{if}0\le x\le 4.5\text{and}5.5\le x\le 10.\end{array}\right.$ (4.3)

$u_3(x,0)=0.$ (4.4)

Figure 2:
Population of wild mosquitoes (a), heterozygous transgenics (b) and homozygous (c) with Dirichlet boundary conditions and initial conditions (4.3) and (4.4).

The amount of mosquitoes obtained at the end of the process results from the integration on the domain Ω of the function obtained in the final time $t=4$. There is no analytic expression for this function, since the curve obtained in $t=4$ comes from a numerical resolution method. The value of this integral was obtained by applying the trapezoid rule, where we obtain $u_1\left(4\right)=44.2654,u_2\left(4\right)=29.5103$ and $u_3\left(4\right)=4.9184$. These values correspond to 56.25% of wild mosquitoes and 43.75% of transgenic (heterozygous and homozygous) mosquitoes, approximately the situation obtained when populations stabilized in Figure 1.

4.2 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.

$u_1(x,0)=u_2(x,0)=\left\{\begin{array}{r}{\left(\sin \left(\pi \left(\frac{x-5.75}{6.75-5.75}\right)\right)\right)}^{100}\text{if}5.75\le x\le 6.75,\\ {\left(\sin \left(\pi \left(\frac{x-4.25}{5.25-4.25}\right)\right)\right)}^{100}\text{if}4.25\le x\le 5.25,\\ 0\text{if}0\le x\le 4.25\text{and}6.75\le x\le 10.\end{array}\right.$ (4.5)

$u_3(x,0)=0.$ (4.6)

Figure 3:
Population of wild mosquitoes (a), heterozygous transgenics (b) and homozygous (c) with Dirichlet boundary conditions and initial conditions (4.5) and (4.6).

4.3 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.

Figure 4 shows the numerical simulation of the dynamics involving wild and transgenic mosquitoes, with initial conditions given by

$u_1\left(0\right)=u_3\left(0\right)=\left\{\begin{array}{r}{\int }_{4.5}^{5.5}10{\left(\sin \left(\pi \left(\frac{x-4.5}{5.5-4.5}\right)\right)\right)}^{100}\text{if}4.5\le x\le 5.5,\\ 0\text{if}0\le x\le 4.5\text{and}5.5\le x\le 10.\end{array}\right.$ (4.7)

$u_2\left(0\right)=0.$ (4.8)

Figure 4:
Population of wild mosquitoes, heterozygous transgenic and homozygous with initial conditions (4.7) and (4.8).

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.

Figure 5:
Population of wild mosquitoes (a), heterozygous transgenic (b) and homozygous (c) with Dirichlet boundary conditions and initial conditions (4.9) and (4.10).

The value obtained for the density of each population at the end of integration, obtained by the trapezoid rule, is $u_1\left(4\right)=19.6735,u_2\left(4\right)=39.3471$ and $u_3\left(4\right)=19.6735$. The proportions obtained are consistent with those obtained in Figure 4.

$u_1(x,0)=u_3(x,0)=\left\{\begin{array}{r}10{\left(\sin \left(\pi \left(\frac{x-4.5}{5.5-4.5}\right)\right)\right)}^{100}\text{if}4.5\le x\le 5.5,\\ 0\text{if}0\le x\le 4,5\text{and}5,5\le x\le 10.\end{array}\right.$ (4.9)

$u_2(x,0)=0.$ (4.10)

5 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.

REFERENCES

• ¹
  F. Catteruccia, T. Nolan, T.G. Loukeris, C. Blass, C. Savakis, F.C. Kafatos & A. Crisanti. Stable germline transformation of the malaria mosquito Anopheles stephensi. Nature, 405 (2000), 959-962.
• ²
  P.R. Couto & S.M. Malta. Interaction between sorption and biodegradation processes in the contaminant transport. Ecological Modelling, 214(1) (2008), 65-73.
• ³
  C. Dufourd & Y. Dumont. Impact of environmental factors on mosquito dispersal in the prospect of sterile insect technique control. Computers & Mathematics with Applications, 66(9) (2013), 1695- 1715.
• ⁴
  R.A. Fisher. The wave of advance of advantageous genes. Annals of eugenics, 7(4) (1937), 355-369.
• ⁵
  V.M. Gantz & E. Bier. The mutagenic chain reaction: a method for converting heterozygous to homozygous mutations. Science, (2015), aaa5945.
• ⁶
  V.M. Gantz , N. Jasinskiene, O. Tatarenkova, A. Fazekas, V.M. Macias, E. Bier & A.A. James. Highly efficient Cas9-mediated gene drive for population modification of the malaria vector mosquito Anopheles stephensi. Proceedings of the National Academy of Sciences, 112(49) (2015), E6736-E6743.
• ⁷
  J. Geiser. “Decomposition methods for differential equations: theory and applications”. CRC Press (2009).
• ⁸
  A.O. H. Diaz, A. A. Ramirez & C. Clavijo. A model for the control of using genetically modified vectors. Journal of Theoretical Biology, 276 (2011), 57-66.
• ⁹
  T.J. Hughes. “The finite element method: linear static and dynamic finite element analysis”. Courier Corporation (2012).
• ¹⁰
  J. Ito, A. Ghosh, L.A. Moreira, E.A. Wimmer & M. Jacobs-Lorena. Transgenic anopheline mosquitoes impaired in transmission of a malaria parasite. Nature , 417 (2002), 452-455.
• ¹¹
  A. Kolmogoroff, I. Petrovsky & N. Piscounoff. Study of the diffusion equation with growth of the quantity of matter and its application to a biology problem. In “Dynamics of Curved Fronts”. Elsevier (1988), pp. 105-130.
• ¹²
  J. Li. Simple mathematical models for interacting wild and transgenic mosquito population. Mathematical Biosciences, 189 (2004), 39-59.
• ¹³
  J. Li. Differential equations models for interacting wild and transgenic mosquito populations. Journal of Biological Dynamics, 2(4-5) (2007), 241-258.
• ¹⁴
  J. Li . Discrete-time models with mosquitoes carrying genetically-modified bacteria. Mathematical Biosciences , 240 (2012), 35-44.
• ¹⁵
  J.S. Lima, A.P. Wyse & A.J.B. Santos. Modelagem da interação entre mosquitos selvagens e transgênicos. In “XIX Encontro Nacional de Modelagem Computacional, XIX EMC e do VII Encontro de Ciência e Tecnologia de Materiais, VII ECTM”. ENMC (2016), p. ENMC.
• ¹⁶
  L.A. Moreira , A.K. Ghosh, E.G. Abraham & M. Jacobs-Lorena. Genetic transformation of mosquitoes: a quest for malaria control. International Journal for Parasitology, 32 (2002), 1599-1605.
• ¹⁷
  L.A. Moreira , J. Wang, F.H. Collins & M. Jacobs-Lorena . Fitness of anopheline mosquitoes expressing transgenes that inhibit Plasmodium development. Genetics, 166(3) (2004), 1337-1341.
• ¹⁸
  C. Noble, B. Adlam, G.M. Church, K.M. Esvelt & M.A. Nowak. Current CRISPR gene drive systems are likely to be highly invasive in wild populations. eLife, 7 (2018), e33423.
• ¹⁹
  T.R. Taha & X. Xu. Parallel Split-Step Fourier Methods for the Coupled Nonlinear Schro¨dinger Type Equations. The Journal of Supercomputing, 32(1) (2005), 5-23. doi:10.1007/s11227-005-0183-5. URL https://doi.org/10.1007/s11227-005-0183-5
  » https://doi.org/10.1007/s11227-005-0183-5» https://doi.org/10.1007/s11227-005-0183-5
• ²⁰
  WHO. World Malaria Report 2014. Technical report (2015).
• ²¹
  WHO. World Malaria Report 2015. Technical report (2016).
• ²²
  WHO. World Malaria Report 2017. Technical report (2017).
• ²³
  A.P. Wyse , L. Bevilacqua & M. Rafikov. Modelagem da interaçao entre mosquitos selvagens e transgênicos. Proceeding Series of the Brazilian Society of Computational and Applied Mathematics, 4(1) (2016).
• ²⁴
  A.P. Wyse , A.J.B. dos Santos, J. dos Santos Azevedo, J.S. de Lima & J.R. de Faria. Modeling the spreading and interaction between wild and transgenic mosquitoes with a random dispersal. PLoS One, 13(10) (2018), e0205879.
• ²⁵
  A.P.P. Wyse. Controle Ótimo do vetor da malária para o modelo matemático sazonal. Laboratório Nacional de Computação Cientifica, (2007).

Publication Dates

• Publication in this collection
  30 Nov 2020
• Date of issue
  Sep-Dec 2020

History

• Received
  10 Dec 2018
• Accepted
  30 Aug 2020