Stability of travelling waves in a Wolbachia invasion

Numerous studies have examined the growth dynamics of Wolbachia within populations and the resultant rate of spatial spread. This spread is typically characterised as a travelling wave with bistable local growth dynamics due to a strong Allee effect generated from cytoplasmic incompatibility. While this rate of spread has been calculated from numerical solutions of reaction-diffusion models, none have examined the spectral stability of such travelling wave solutions. In this study we analyse the stability of a travelling wave solution generated by the reaction-diffusion model of Chan&Kim (2013) by computing the essential and point spectrum of the linearised operator arising in the model. The point spectrum is computed via an Evans function using the compound matrix method, whereby we find that it has no roots with positive real part. Moreover, the essential spectrum lies strictly in the left half plane. Thus, we find that the travelling wave solution found by Chan&Kim (2013) corresponding to competition between Wolbachia-infected and -uninfected mosquitoes is linearly stable. We employ a dimension counting argument to suggest that, under realistic conditions, the wavespeed corresponding to such a solution is unique.


Introduction
Wolbachia are common endosymbiotic bacteria that are estimated to infect up to 66% of all insect species [10]. They are primarily vertically transmitted and can induce reproductive phenotypes in its host to confer a reproductive advantage and hence improve chances of persistence. A well studied induced reproductive phenotype is cytoplasmic incompatibility (CI), whereby offspring of infected males and uninfected females have an increased probability of death. This has lead to the proposal of a deliberate Wolbachia introduction into wild mosquito populations to reduce transmission of vector-borne diseases, since particular CI-inducing Wolbachia strains have been shown to reduce proliferation of various viruses (see Brelsfoard & Dobson [3] and references therein). Two particular CI-inducing strains of Wolbachia, WMel and WMelPop, have received much attention due to evidence of these strains inhibiting dengue transmission in Aedes Aegypti mosquitoes.
Although CI has been proposed as a key mechanism for the success of Wolbachia-based strategies, fitness reducing phenotypes are also a result of certain Wolbachia strains. For example, both WMel and WMelPop reduce both fecundity and lifespan of infected females [15,22,23]. Turelli [20] first showed that the amalgamation of these Wolbachia-induced effects result in a strong Allee effect for an invasion, that is, infection is only successful if initial infection densities are above the critical Allee threshold. This threshold has been found by various modelling studies to be dependent on the fecundity cost of infection, the strength of CI in reducing offspring from infected males and uninfected females, and the probability of successful Wolbachia transmission [4,8,12,16,21].
Many studies have analysed the growth dynamics of a deliberate Wolbachia introduction via non-spatial models, but relatively few studies have studied the dynamics of the spatial spread (see Hancock & Godfray [9] and Barton & Turelli [2] for examples). This study is motivated by the reaction-diffusion model of Chan & Kim [4], who examine the spatial spread of Wolbachia in a homogeneous environment by incorporating slow and fast compartments in their model. Chan & Kim [4] numerically show that there exists a travelling wave solution corresponding to a Wolbachia invasion, that is, competition between Wolbachia-infected and -uninfected mosquitoes, and estimate the wavespeed corresponding to this solution. A simplified model which uses a weighted average of the slow and fast diffusion coefficients was found to yield similar wavespeeds. This model assumes perfect vertical transmission of Wolbachia, an increased fecundity for infected individuals, longevity reducing effects of 10% and that Wolbachia infection induces CI. Here we perform a spectral stability analysis by examining the linearised operator arising in this simplified model. We show that the essential spectrum is bounded to the left-half plane for all relevant biological parameter values and show that the point spectrum contains no elements in the right-half plane. Moreover, we show that there exists a travelling wave solution for only a unique wavespeed.

Problem setup
The non-dimensionalised system of partial differential equations from Chan & Kim [4] is given by where u is the density of Wolbachia infected mosquitoes, v is the density of uninfected mosquitoes, S = u + v and A = u/S. Following Chan & Kim [4], we let s h = 0.45, F = 1.0526 and µ = 0.0162, which correspond to parameter values for Aedes aegypti at 30 • C. The parameters F −1 , µ, ρ and s h correspond to the relative fecundity of infected females to uninfected females, mortality rate, advection rate and probability of embryo death due to cytoplasmic incompatibility, respectively. Additionally, we let α = 1.1, which reflects the 10% relative reduction in lifespan associated with the WMel strain of Wolbachia.
We linearise about the travelling wave solution (û(z),v(z)) via the substitution where p(z, t) and q(z, t) are perturbations in H 1 , ∀t ∈ R. Collecting first order perturbation terms, we obtain We define the linear operator L by which has the corresponding eigenvalue problem (L − λ) p q = 0. We introduce the substitutions s = p z and t = q z to convert the eigenvalue problem into a first order boundary value problem, and denote this equivalent operator of (L − λ) p q = 0 by T (p, q, s, t) T , where T (y)(z) = d dz − A(z, λ) y and y = (p, q, s, t) T . This process yields and To assess the stability of travelling wave solutionû andv we need to locate the spectrum of the linearised operator L as an operator on H 1 × H 1 . If (L − λ) −1 does not exist or is unbounded for λ ∈ C, then λ is in the (ii) The point spectrum, defined by We define A ± (λ) := lim z→±∞ A(z, λ), which are given by The asymptotic operator of T is given by where A ∞ is the piecewise spatially constant matrix (2.11)

Essential spectrum
The PDE given in (2.2) is autonomous, and so the only non-constant coefficients are due to the functionsû(z) andv(z) in the reaction terms. These are heteroclinic orbits in phase space (connecting e − = (1 − αµ, 0) and e + = (0, 1 − µ F )) which decay exponentially as z → ±∞ as shown in Section 5. This shows that L is exponentially asymptotic. From Kapitula & Promislow [11, Theorem 3.1.11] it follows that L is a relatively compact perturbation of the asymptotic operator L ∞ , defined as the limit of L as z → ±∞ and equivalent to the operator T ∞ (0). Then by Weyl's Essential Spectrum Theorem, we have that σ ess (L) = σ ess (L ∞ ), or A crucial concept behind the spectrum of an operator is the existence of an exponential dichotomy. Essentially, this states that each solution to (2.6) decays exponentially either in forward or backward z. For spatially constant matrices, the existence of an exponential dichotomy simply means that the matrix is hyperbolic. We define the Morse index of a constant matrix A to be the dimension of the unstable subspace associated with A, and let i ± (λ) denote the Morse indices of the asymptotic matrices A ± , given by Eq. (2.8) and Eq. (2.9). It can be shown that for λ ∈ C such that T ∞ is Fredholm, we have ind( Thus, we can characterise the essential spectrum of L ∞ as where E c denotes the center subspace associated with the asymptotic linearised system. The spatial eigenvalues of A − (λ) and A + (λ) are respectively given by These spatial eigenvalues are non-hyperbolic when η ± = ik. Substituting this into Eq. (3.2) and (3.3) respectively yields the dispersion relations These form four parabolas in the complex plane parametrised by k. For λ in between the region bounded by has three stable eigenvalues and one unstable eigenvalue; to the left of the region A − (λ) has four stable eigenvalues and to the right of the region A − (λ) has two stable and two unstable eigenvalues. This is also true for λ 1,2 + (k) and A + (λ). Thus the essential spectrum is given by the region bounded between λ 1 − and λ 1 + , and also between λ 2 − and λ 2 + ; this is shown in Figure 1.

Re λ
Im λ Figure 1: The essential spectrum of L is given by λ in the shaded regions. The blue dashed and solid lines represent λ 1,2 − and λ 1,2 + respectively. The red line indicates the absolute spectrum given by Eq. (3.6) and the red dot at the origin represents an eigenvalue. Note that these are not drawn to scale for visualisation purposes.
It will be convenient for us later on to know the location of the so-called absolute spectrum. The absolute spectrum is not spectrum per se, but its location characterises the breakdown of the analytic continuation (in terms of the spectral parameter λ) of the stable and unstable eigenspaces of the matrices A ± (λ). It thus follows that the absolute spectrum coincides with a branch cut of the Evans function [11,Section 3.2]. For the case at hand, the absolute spectrum can be defined [19] as the set in the complex plane where a pair of the eigenvalues of A + (λ) have equal real parts. This is the set We note that from Eq. (3.6), (3.4) and (3.5), the continuous and absolute spectrum are always bounded to the left-half plane for biologically relevant parameter constraints F > 1, µ > 0, s h > 0 and α > 1.

Point spectrum
The existence of a travelling wave solution implies a heteroclinic connection in (2.4) between the equilibria e − = (1 − αµ, 0) and e + = (0, 1 − µ F ). We denote the unstable subspace of the matrix A − by U − and the stable subspace of the matrix A + by S + . To the right of the essential spectrum, we have that the dimension of U − , which we denote by k, and dimension of S + sum to 4, the dimension of the entire phase space. For our case, k = 2.
We initialise Eq. (2.6) at z = −∞ with ζ − 1 , ζ − 2 and at z = ∞ with ζ + 1 , ζ + 2 , and solve the system towards a matching point, which we pick to be z = 0. We denote the solutions of the former by w − 1 (0, λ), w − 2 (0, λ) and the latter by w + The Evans function is defined by

Compound matrix method
The Evans function is numerically difficult to compute due to the stiffness of the problem stemming from the difficulty of resolving different modes of growth and decay. For example, since we have that η − 1 > η − 2 , any numerical errors occurring when solving for w − 2 (0, λ) will grow at a rate proportionate to exp(η − 1 z). Thus, although the solutions w − 1 (z, λ), w − 2 (z, λ) are linearly independent at z = −∞, they quickly become numerically linearly dependent. Several methods have been proposed to overcome this issue such as the method of continuous orthogonalisation, the compound matrix method, Magnus methods and Grassmanian spectral shooting [6,7,13,17,18]. Following Allen & Bridges [1], we employ the compound matrix method which converts the problem into the six-dimensional wedge product space ∧ 2 (C 4 ), with basis B = {e 1 ∧e 2 , e 1 ∧e 3 , e 1 ∧e 4 , e 2 ∧e 3 , e 2 ∧e 4 , e 3 ∧e 4 }, where {e 1 , e 2 , e 3 , e 4 } is the standard basis for C 4 . The numerical advantage of this approach is that the evolution of w − 1 , w − 2 and w + 1 , w + 2 are incorporated into a single trajectory given by w − 1 ∧ w − 2 and w + 1 ∧ w + 2 respectively. The coordinate vector of w − 1 ∧ w − 2 relative to the basis B is given by and the second subscript in w − i,j denotes the jth element within the vector w It can be shown (see Allen & Bridges [1]) that φ(z) = φ − (z), φ + (z) satisfy the equation whereÃ is the induced matrix given bỹ with A ij given by Eq. (2.7).
We initialise the problem at z = −∞ with φ − (−∞) and at z = ∞ with φ + (∞) and solve towards the matching point z = 0. The Evans function is then defined to be (4.9) For numerical stability, we scale the solution according to its exponential growth and decay rates by letting (4.11) where the Evans function is equivalent to Eq. (4.9) except with φ replaced with ψ.
Since the Evans function is analytic to the right of the essential spectrum, we have via the Argument Principle (4.12) where N is the number of zeroes in the interior of the region enclosed by C.
To check for eigenvalues of L with positive real part, we set up a closed semi-circle contour C excluding the origin, as shown below in Figure 2. We let r s and r b denote the radius of the smaller and larger circular arc respectively.

Re λ
Im λ     respectively. The only roots of the Evans function are at λ = 0 and λ ≈ −0.0026075, the latter being the edge of the absolute spectrum, which we denote by γ A . The former is due to translational invariance of the travelling wave solution (û,v), while the latter is due to γ A being a branch point of D(λ). We note that the root at λ = 0 is simple, and because of translational invariance must persist throughout all nearby parameter regimes. We have already shown that the root of the Evans function at the branch point is in the left half plane and thus we conclude that no new eigenvalues can be introduced by perturbation.
Since system (2.1) has no spectrum in the right half plane, the solution given byû andv is spectrally stable. Moreover, as the linearised operator L is an exponentially asymptotic operator, we have that it is also a sectorial operator [5,11, see Chapter XVII, §6, Proposition 3 in the former and Example 4.1.8 in the latter].
Thus, spectral stability of the travelling wave solutionû andv also implies linear stability. We refer the reader to Section 5 for details on computingû andv.  Figure 6: Plot of the Evans function given by Eq. (4.9). The blue and red solid lines in both plots show D(λ) with α = 1.1 and α = 1 respectively. The dashed lines in plot (a) mark the edge of the absolute spectrum corresponding to each α. For α = 1.1, the only roots are at λ = 0 and at the edge of the absolute spectrum λ = −0.002607, whereas for α = 1 we were unable to detect a zero at the edge of the absolute spectrum due to its proximity to the origin.
We note that for parameter values such that γ A lies closer to the origin, for example α = 1 (corresponding to γ A ≈ −0.001435), the method described above fails to detect a zero for the Evans function evaluated at γ A , although the qualitative behaviour in the right-half plane remains the same. We show this in Figure 6.

Wave profile
To compute A(z, λ) explicitly at any z requires either a numerical or exact solution forû(z) andv(z) satisfying Eq. (2.1). We use MATLAB's bvp4c solver to find a numerical solution corresponding to the case where the populations represented byû andv are in competition. The boundary conditions listed in (2.6) are not sufficient for bvp4c to find a unique solution. We note that since the linearisation of (2.2) as z → ±∞ is given by y = A ∞ (0)y, we have with η ± 2 as defined in Eq. (4.1)-(4.2), but with λ = 0. To ensure uniqueness of the solution, we include this information on the derivatives in the boundary conditions by setting where L is a large number chosen to represent numerical infinity (chosen to be L = 200 in Figure 7). The numerical solution of the wave profilesû andv are shown in Figure 7. Due to the bistability in the system, and from the spatial dynamics of the PDE u t = u xx + u(1 − u)(u − a) (essentially a reaction-diffusion equation with bistable reaction term/strong Allee growth dynamics, see Lewis & Kareiva [14]), we expect that there is a unique wavespeed c * , which we numerically determine to be approximately 0.027, for which there is a heteroclinic connection between e − and e + in system (2.2).
Linearising about e − and e + shows that the dimensions of W u (e − ) (the unstable manifold of e − ) and W s (e + ) (the stable manifold of e + ) are both equal to 2. We assume that these two manifolds are generic and intersect tranversely, which we interpret as a biologically realistic assumption. Then we consider extending system (2.2) with c = 0. Due to transversality, we have that codim(W u (e − ) ∩ W s (e + )) = codim(W u (e − )) + codim(W s (e + )) = 2 + 2, which leads to dim(W u (e − ) ∩ W s (e + )) = 1. Figure 8 provides a schematic of this argument.

Discussion
In this study we have shown the linear stability of a travelling wave solution to a model for Wolbachia spread. This is achieved by computing the essential, point and absolute spectrum of the linearised operator and showing the absence of spectrum in the right-half plane. We prove that the essential and absolute spectrum is bounded to the left-half plane for all biologically relevant parameter settings. Due to the numerical nature of locating the point spectrum, we only show that there is no point spectrum in the right-half plane for fixed parameter values.
Our results suggest that although Wolbachia may be difficult to establish in a local area due to a CI-induced strong Allee effect in the growth dynamics, once it is established the spread of infection is a stable phenomenon.
In addition to our study being an investigation of Wolbachia spread dynamics, we present our study as an example of a dynamical systems approach to determining stability of travelling wave solutions in a system of PDEs. Models demonstrating such solutions are ubiquitous, particularly as mathematical modelling is becoming increasingly integrated with the scientific method. The dynamical systems tools we have used in this study can be applied to a wide variety of models currently used in mathematical biology; we believe that their application will improve understanding of the dynamics generated by a model and motivate research in more complicated biological models.
One of the key obstacles impeding the wider use of the tools in this study is the difficulty in computing the point spectrum via the Evans function. There are two key difficulties regarding this. Firstly, numerical methods for evaluating the Evans function sometimes fail, due to the evaluation requiring the solution to a stiff problem.
Although in this study we have successfully used the compound matrix method, it is not guaranteed to work for all cases. Secondly, evaluating the Evans function requires the solution whose stability we are interested in. While obtaining such a solution is straightforward when the system of PDEs is exactly solvable, in many cases it is not exactly solvable and instead one must rely on a numerical solution obtained through solving a boundary value problem (see Section 5). Depending on the dimensionality of the problem and the dimensions of the stable and unstable subspaces at the equilibria, this can be a non-trivial numerical problem.