A wolbachia infection model with free boundary

Scientists have been seeking ways to use Wolbachia to eliminate the mosquitoes that spread human diseases. Could Wolbachia be the determining factor in controlling the mosquito-borne infectious diseases? To answer this question mathematically, we develop a reaction-diffusion model with free boundary in a one-dimensional environment. We divide the female mosquito population into two groups: one is the uninfected mosquito population that grows in the whole region while the other is the mosquito population infected with Wolbachia that occupies a finite small region and invades the environment with a spreading front governed by a free boundary satisfying the well-known one-phase Stefan condition. For the resulting free boundary problem, we establish criteria under which spreading and vanishing occur. Our results provide useful insights on designing a feasible mosquito releasing strategy to invade the whole mosquito population with Wolbachia infection and thus eventually eradicate the mosquito-borne diseases.


Introduction
Recently, several public health projects have been launched, in China [28], USA [25] and France [22], with an aim to fight mosquito populations that transmit Zika virus, Dengue fever and Chikungunya. All of these projects involve the release of male Aedes aegypti mosquitoes infected with the Wolbachia bacteria to the wild. For instance, 20000 male Aedes aegypti mosquitoes carrying Wolbachia bacteria were released on Stock Island of the Florida Keys in the week of April 20, 2017. Google's Verily is about to release 20 million machinereared Wolbachia-infected mosquitoes in Fresno (see [25]). A factory in Southern China is manufacturing millions of "mosquito warriors" (male Aedes aegypti mosquitoes carrying Wolbachia bacteria) to combat epidemics transmitted by mosquitoes [28].
The science behind these projects is based on the following two facts: (i) Wolbachia often induces cytoplasmic incompatibility (CI) which leads to early embryonic death when Wolbachia-infected males mate with uninfected females and (ii) Wolbachia-infected females produce viable embryos after mating with either infected or uninfected males, resulting in a reproductive advantage over uninfected females. In practice, Wolbachia has been successfully transferred into Aedes aegypti or Aedes albopictus by embryonic microinjections, and the injected infection has been stably maintained with complete CI and nearly perfect maternal transmission [1,12,16,17,23,33,34,36]. Thus, the bacterium is expected to invade host population easily driving the host population to decline. Successful Wolbachia invasion in Aedes aegypti has been observed by Xi et al. in the laboratory caged population within seven generations [35].
By releasing Aedes albopictus mosquitoes infected with Wolbachia bacteria into the wild, it is expected that over a long time period, the wild Aedes aegypti mosquito population would decline drastically and hopefully be completely replaced by infected mosquitoes so that the mosquito-borne infectious diseases such as Zika, Dengue fever and Chikungunya would be eradicated. To qualitatively examine if Wolbachia can effectively invade the wild uninfected mosquito population, Zheng, Tang and Yu [40] considered the following model: where u denotes the number of reproductive infected insects and v denotes uninfected ones, b 1 and b 2 denote half of the constant birth rates for the infected and uninfected insects respectively. The parameter δ 1 (resp. δ 2 ) denotes the density-dependent death rate for the infected (resp. uninfected) population. The birth rate of uninfected mosquitoes is diminished by the factor v u+v due to the sterility caused by cytoplasmic incompatibility (CI) for mating between infected males and uninfected females.
Let us now recall the origin of system (1) with some details. Let r f and r m denote the number of released female mosquitoes and the number of released males respectively and suppose they were infected with Wolbachia. Also, assume that r f and r m satisfy where T (t) = r f + r m + I f + I m + U f + U m denotes the total population size, with U f , U m , I f and I m standing for the numbers of uninfected reproductive females, uninfected reproductive males, and infected reproductive females and males other than those from releasing, respectively. Let b I (resp. b U ) be the natural birth rate of the infected (resp. uninfected) mosquitos and 0 ≤ δ ≤ 1 be the proportion of mosquitos born female. Then the proportion of mosquitos born male is 1 − δ. With complete CI (see Table 1) and perfect maternal transmission, we have One can easily verify that both r f and r m approach 0 as t → +∞. We denote by Assuming equal determination case, which means that δ = 1/2, I f = I m and U f = U m , then system (1) can be obtained by setting b 1 = b I /2 and b 2 = b U /2. In order to obtain the spatiotemporal dynamics of (1), Huang et al. [13,14] studied the following reaction-diffusion system: In (5), d 1 and d 2 are the diffusion rates, ∆ denotes the Laplace operator in the spatial variable x, and ν denotes the unit outward normal vector to the boundary of Ω. We mention that (5) is obtained from a delay differential equation model in [40] after ignoring the delay factor and incorporating the spatial inhomogeneity.
Similarly, there has been several mathematical models formulated to describe the Wolbachia spreading dynamics [15,37,38,41]. These models focused on studying the subtle relation between the threshold releasing level for Wolbachiainfected mosquitoes and several important parameters including the CI intensity and the fecundity cost of Wolbachia infection. We also note that female Aedes aegypti mosquitoes infected with the Wolbachia bacteria were initially released at a specific site. Hence, the infected female mosquitoes initially occupy only a small region, while the wild uninfected females are distributed over the whole area.
To model the spatial spreading of Wolbachia in the wild mosquito population and explore the possibility that the infection can indeed occupy the whole region, it is natural to consider system (5) under the setting of a free boundary problem.
In this work, we consider the following free boundary problem in one-dimensional space: The equation governing the movement of the spreading front x = h(t) is deduced in a manner similar to that in Section 1.3 of [2]. It is known as the one-phase Stefan condition in the literature. This type of free boundary condition has been widely used in previous works such as [4,5,6,9,18,19,20,24,29,31,32]. We will first analyze system (6) with constant birth rates b 1 and b 2 in Section 3. Note that environmental variables such as especially available water surfaces and humidity have huge impacts on birth rates [7], we also extend our study to the case with space-dependent birth rates b 1 (x) and b 2 (x) in Section 4, while the natural death rate is assumed to be spatially independent for simplicity.
Throughout this paper, we assume that b 1 (x) and b 2 (x) satisfy the following conditions, unless otherwise stated: C 0,θ ([0, +∞)) is the Hölder space with Hölder exponent θ. The initial conditions u 0 and v 0 are assumed to be bounded and satisfy For the free boundary problem (6)-(7), the main question we are concerned about is whether the infected population can eventually occupy the whole space or not. The main goal of this work is to derive conditions under which the spreading occurs. If spreading occurs, then the whole mosquito population will become infected with Wolbachia bacteria and this leads to the extinction of the mosquito population and eventually the eradication of mosquito-borne diseases.
Organization of the paper. The paper is organized as follows. We first establish the global existence and uniqueness of solutions to the free boundary problem (6) in Section 2. In Section 3, we present a detailed analysis of a specific case of model (6). In Section 4, we study the population dynamics of infected mosquitoes in a heterogeneous environment with a free boundary condition. In order to better understand the effects of dispersal and spatial variations on the outcome of the competition, we study system (6) over a bounded domain with Neumann boundary conditions. We summarize our results in the last section.
The next result provides some bounds on the solutions to system (6) with initial conditions (7).
Thus, h ′ (t) > 0 for t ∈ (0, T ]. Next, we consider the initial value problem From the comparison principle, we know that Similarly, we can show that To prove (iii), we first consider the auxiliary function , We also note that and Applying the comparison principle, we get Since ω 1 (t, h(t)) = 0 = u(t, h(t)), we then have Bearing the above result in mind, we can show that the local solution obtained in Theorem 2.1 can indeed be extended to all t > 0. Proof. Let [0, T max ) be the maximal time interval in which the unique solution exists. We will show that T max = ∞. Suppose to the contrary that T max < ∞. In view of Lemma 2.2, there exists positive constants M 1 , M 2 and Λ, independent of T max , such that for t ∈ [0, T max ], Fix δ ∈ (0, T max ) and K > T max . Using the standard L p estimates together with the Sobolev embedding theorem and the Hölder estimates for parabolic equations (see Lunardi [21] for eg.), we can find M 3 depending only on δ, K, M 1 and M 2 such that where we used the convention that u(t, x) = 0 for x ≥ h(t). By virtue of the proof of Theorem 2.1 in [10], there exists a τ > 0 depending only on M 1 , M 2 and M 3 such that the solution of (6) with the initial time T max − τ 2 can be extended uniquely to the time T max + τ 2 , which contradicts the definition of T max . Thus, T max = +∞ and the proof is complete.
3 The special case of constant birth rates System (5) was investigated in [13,14] for two disjoint cases. Namely, the fitness benefit case and the fitness cost case. Define κ 1 and κ 2 as Wolbachia is said to have the fitness benefit if κ 1 > κ 2 , which means that the local area is more (or at least equally) favourable for infected mosquitoes. The fitness cost case is represented by κ 1 < κ 2 , see [40].
In this section, we assume that b i (x) = b i for i = 1, 2, where b i are positive constants. In other words, we have the constant-coefficient free boundary problem given by System (11) is essentially a competition model. For the fitness benefit case, κ 1 > κ 2 , u is the so-called superior competitor and v the inferior competitor (see [10]). For the fitness cost case, κ 1 < κ 2 , (11) represents a strong competition [26]. Throughout this section, we always assume u is a superior competitor. That is, the Wolbachia infection has a fitness benefit. The strong competition case is usually more complicated to be studied mathematically. To the best of our knowledge, results for competition models with a free boundary are very limited in strong competition case. Further details can be seen in [42,43]. We organize this section as follows. In subsection 3.1 we present some preliminary results, which play a role in proving our main results. Subsection 3.2 is devoted to the vanishing case. The invasion dynamics is studied in detail in Subsection 3.3. A rough estimation of asymptotic spreading speed of Wolbachia invasion is given in Subsection 3.4.

Preliminary results
The following result holds.
We recall the following comparison principle.
Lemma 3.2 (Comparison principle [10]). Assume that 0 ≤ T 0 < T < +∞ and h, h ∈ C 1 ([T 0 , T ]). Denote by and Let (u, v, h) be the unique solution of (11). Then, The following follows from Lemmas A.2 and A.3 in [39]. with (b) Let a, b and q be fixed positive constants. For any given ε > 0 and L > 0, where V (t, x) is a continuous and non-negative function satisfying and We are now in the position to present part of our main results.

The vanishing case
We consider the vanishing case in this subsection.
be the solution of system (11) with initial data (7).
The criteria for spreading and vanishing are given in the following theorem.
Proof. Note that h(t) is nondecreasing. We only need to show that h ∞ < +∞ implies h ∞ ≤ h * 0 . It follows from Theorem 3.4 that h ∞ < +∞ implies Let u be the solution of the following problem By the comparison principle, we have u(t, x) ≤ u(t, x) for all t ≥ T and x ∈ which is a contradiction. Therefore, h ∞ ≤ h * 0 and this completes the proof.
It follows from the comparison principle that Clearly, u(T 1 , x) is independent of µ. Now, we consider the following system.
By Theorems 3.6 and 3.7, we can also derive spreading criteria in terms of the diffusion coefficient d 1 , for any fixed h 0 .
where h 0 is any prefixed positive constant. Then, spreading occurs provided that either Our next result is a criterion on "vanishing".
Proof. Consider the following problem Lemma 3.2 applies and yields that h(t) ≤h(t) and u(t, x) ≤ū(t, x) for t > 0 and 0 ≤ x ≤ h(t).

The spreading speed
If spreading occurs, it is important to estimate the spreading speed of h(t).
Following an idea in [11], one can obtain a rough estimate of the spreading speed as stated in the following theorem.
where s * is the minimal speed of the traveling waves to the problem related with (11) in the entire space. This estimation of the spreading speed is independent of µ.
However, in the fitness benefit case, we can derive an estimate better than the one in Theorem 3.4. We first recall Proposition 5.1 of [10].
Proof. Note that Thus, the pair (u, h) is a subsolution to the problem By the comparison principle, h(t) ≤h(t) for t > 0. Theorem 4.2 of [9] yields that Note that lim sup t→+∞ v(t, x) ≤ κ 2 uniformly for x ∈ [0, +∞) and h ∞ = +∞. Then, Next, we consider the following problem By the comparison principle, we obtain h(t) ≥ h(t) for t > T ε . From Theorem 3.6, we know that h(∞) = +∞. Using a similar argument as above, we have

The free boundary problem with a heterogeneous birth rate
In this section, we consider the free boundary problem (6)-(7) with the heterogeneous birth rates b 1 (x) and b 2 (x).

Some useful lemmas
In this subsection, we first study a related eigenvalue problem: Problem (37) admits a positive principal eigenvalue λ 1 determined by We state two hypotheses that we refer to when needed. We use a generic symbol B(x) in the statement of the hypotheses. The function B(x) will be replaced accordingly (by b, b 1 or b 2 ) in the rest of this Section.
Remark 1. In order to compare the principal eigenvalues λ 1 associated with different parameters, we denote the principal eigenvalue λ 1 by λ 1 (d, h 0 ). When we fix h 0 and study the property of λ 1 as d varies, we write λ 1 = λ 1 (d). Similarly, we write λ 1 = λ 1 (h 0 ) when d is fixed while h 0 varies.
We gather the following known results about the dependance of λ 1 on d and h. (i) λ 1 (d) is increasing with respect to d.  replaced by b(x). Then, λ 1 = λ 1 (h 0 ) has the following properties: For the reader's convenience, we also recall some facts related to the following problem The proof of the next lemma follows from Lemma 5.2 and Lemma 6.2 of [44]. where φ v * is the unique positive solution of the following elliptic problem

Sharp criteria for spreading and vanishing
Let us first consider the vanishing case. The proof is similar to that of Theorem 3.4, above.
In order to obtain sharp criteria for spreading, we require stronger conditions on b 1 (x) and δ 1 . Namely, we assume that Our assumption (41) is not excessive in the sense that, when b i and δ i (i = 1, 2) are constant, we have φ v * = b2 δ2 . Consequently, b 1 − δ 1 φ v * is a positive constant over the interval [0, +∞).
Proof. First, we consider the following equation: Since b 2 (x) satisfies the hypotheses of Lemma 4.3, all solutions of (42) with non-trivial non-negative initial values converge to φ v * as t → ∞. It follows, from the comparison principle, that v ≤ v for all t > 0 and x > 0. Since lim t→+∞ v(t, x) = φ v * uniformly in any compact subset of [0, ∞), then, for any ε > 0, there exists T > 0 such that v(t, x) ≤ φ v * + ε, for t ≥ T . Consider the following eigenvalue problem: It is well known that the principal eigenvalue λ 1 can be characterized by Using (iii) of Lemma 4.1, for any fixed h 0 , there exists d * 1 such that In this theorem, we have 0 < d 1 < d * 1 . Let us set u = δϕ 1 (x), for t ≥ T and x ∈ [0, h 0 ] (here ϕ 1 (x) is the corresponding eigenfunction of λ 1 ). Choose δ > 0, small enough, so that A straightforward calculation leads to  Proof. Similarly, we consider the following equation Since b 2 (x) satisfies the hypotheses of Lemma 4.3, all solutions of (45) with nontrivial and nonnegative initial conditions converge to φ v * as t → ∞. It follows from the comparison principle that v ≤ v for t > 0, x > 0. Since lim sup t→+∞ v(t, x) = φ v * uniformly in any compact subset of [0, ∞). So for any Consider the following eigenvalue problem: The principal eigenvalue λ 1 is characterized by is the corresponding eigenfunction of λ 1 ). Choose δ > 0 small enough so that By the comparison principle, Similarly, we have h ∞ = +∞; hence, spreading occurs.
Theorem 4.7. If d 1 > d * 1 and u 0 is small enough, then "vanishing" occurs. Proof. We consider the following problem as an auxiliary to the first equation of (6): Denote the principal eigenvalue λ 1 and the corresponding positive eigenfunction One can verify that there exists d * 1 such that λ 1 > 0, when d 1 > d * 1 . Furthermore, it follows, from Theorem 4.2 in [44], that there exists a constant B such that φ ′ 1 (x) ≤ 2h 0 Bφ 1 (x) for all x ∈ [0, h 0 ]. Now, we can use the following auxiliary functions, which were constructed in [44]. Let , for t ≥ 0 and The conditions on α and β will be determined later. If we let 0 < α ≤ 1, direct calculations show that Since h(t) → h 0 as α → 0, we can find sufficiently small α 1 , such that Moreover, there exists α 2 > 0, small enough, such that , for α ≤ α 2 .
Let α = min 1, λ1 4 , α 1 , α 2 . Direct calculation leads to . Then, In order to apply the comparison principle, we choose u 0 small enough such that Form the comparison principle, we have h(t) ≤ h(t) for t > 0 and u(t, x) ≤ u(x, t) for t > 0 and x ∈ [0, h(t)].
Moreover, we can derive vanishing criteria in terms of the coefficient µ when d 1 > d * 1 .

Thus,
This implies that vanishing occurs.
Next, we will prove the following conclusions.
It follows from the comparison principle that 0 ≤ u ≤ū for t > 0 and x ∈ (0, l).
Since l ≤ π 2 Under some assumptions, stated below, we can obtain the asymptotic spreading speed from Theorem 3.6 of [8].

Summary and conclusions
We studied a reaction-diffusion model with a free boundary in one-dimensional environment. The model is developed to better understand the dynamics of Wolbachia infection under the assumptions supported by recent experiments such as perfect maternal transmission and complete CI.
In the special case of constant birth rates, we only considered the fitness benefit case. For the fitness benefit case, where the environment is more favorable for infected mosquitoes, our results show that the spreading of Wolbachia infection occurs if either the size of the initial habitat of infected population h 0 is large enough, say h 0 ≥ h * 0 (Theorem 3.6), or the boundary moving coefficient µ is sufficiently large (µ ≥μ) in case of h 0 < h * 0 (Theorem 3.7). A rough estimate on the spreading speed of h(t) is also provided. Moreover, if h 0 < π 2 d1 b1 < h * 0 and µ ≤ µ, then the infection cannot spread and h ∞ < +∞.
The case of inhomogeneous (spatially dependent) birth rates is treated in Section 4. Detailed criteria for spreading and vanishing are derived in Subsection 4.2 with the aid of spectral properties of relevant eigenvalue problems.