*Article* **Individual Variability in Dispersal and Invasion Speed**

#### **Aled Morris 1,2, Luca Börger 1,3 and Elaine Crooks 1,2,\***


Received: 31 July 2019; Accepted: 22 August 2019; Published: 1 September 2019

**Abstract:** We model the growth, dispersal and mutation of two phenotypes of a species using reaction–diffusion equations, focusing on the biologically realistic case of small mutation rates. Having verified that the addition of a small linear mutation term to a Lotka–Volterra system limits it to only two steady states in the case of weak competition, an unstable extinction state and a stable coexistence state, we exploit the fact that the spreading speed of the system is known to be linearly determinate to show that the spreading speed is a nonincreasing function of the mutation rate, so that greater mixing between phenotypes leads to slower propagation. We also find the ratio at which the phenotypes occur at the leading edge in the limit of vanishing mutation.

**Keywords:** invasive species; linear determinacy; population growth; mutation; spreading speeds; travelling waves

#### **1. Introduction**

The speed at which a species expands its range is a fundamental parameter in ecology, evolution and conservation biology. Knowledge of this speed enables us to predict the ability of a species to keep up with the rate at which the climate changes or the rate at which an exotic species invades, representing two prominent ecological challenges [1,2]. It is known that traits such as dispersal and population growth affect the rate at which a species expands its range, and there has been a suggestion in recent work that polymorphism in traits could cause a species invasion to occur at a faster rate than a single morph would in isolation [3,4]. Understanding the effect that each trait of a species has, and could potentially have, on its rate of spread is therefore important to understanding how the spread of a species can evolve.

Most common models of invasions in population dynamics incorporate aspects of dispersal and growth, e.g., works by the authors of [5–8], however the mutation of one phenotype to another has been a less common inclusion. Even the addition of a simple mutation term can dramatically affect the behaviour of a model. A review of Cosner [9] singles out two models that involve mutation and multiple dispersal strategies in a population of a species: the model introduced in Elliott and Cornell [3] to investigate dispersal polymorphism for two morphs, in which a simple linear mutation is used, and the model of Bouin et al. [10], motivated by the destructive invasion of cane toads across northern Australia, in which mutations are considered to act as a diffusion process in the phenotype space.

Elliott and Cornell assume that the spread rate of the two phenotypes in their system, usually referred to as the spreading speed, is determined by the linearisation of their system at the extinction state zero. This assumption can be rigorously proved to hold under reasonable conditions on the parameters; see the framework of Girardin [11], which applies in particular to this model, and also, for an alternative approach, Morris [12], where the assumption is proved using earlier results of Wang [13] in the case where the mutation rate is small, which is generally the case for all organisms since natural selection typically

acts to minimize mutation rate [14]. Moreover, in addition to being linearly determined, the spreading speed equals the minimal speed of a class of travelling waves [11,12], mimicking well-known results on travelling waves and spreading speeds for the Fisher–KPP equation [15–17].

Knowing that the rate of spread is linearly determinate and linked to travelling wave speeds provides a powerful tool that we exploit here to deduce ecologically-important information about the invasion of trait-structured species using the model introduced in [3]. In particular, we establish results on the dependence of spreading speeds on the mutation rate, and on the composition of the leading edge of minimal speed travelling waves in the limit of vanishing mutation. Some of our results focus, as in [3], on the case when different morphs have varying dispersal abilities or strategies, and in addition, there is a trade-off between dispersal and growth. Such trade-offs are exhibited by many species, including certain plants, insects and terrestrial arthropods; see, for instance, the review of Bonte et al. [18] on the costs of dispersal.

Elliott and Cornell's model examines the interaction between an establisher phenotype with population density *ne*, and a disperser phenotype with population density *nd*, using a Lotka–Volterra competition system:

$$\begin{aligned} \frac{\partial n\_{\varepsilon}}{\partial t} &= D\_{\varepsilon} \frac{\partial^2 n\_{\varepsilon}}{\partial \mathbf{x}^2} + r\_{\varepsilon} n\_{\varepsilon} (1 - m\_{\text{cc}} n\_{\varepsilon} - m\_{\text{cd}} n\_d) - \mu \epsilon n\_{\varepsilon} + \mu d n\_d\\ \frac{\partial n\_d}{\partial t} &= D\_d \frac{\partial^2 n\_d}{\partial \mathbf{x}^2} + r\_d n\_d (1 - m\_{\text{dc}} n\_{\text{c}} - m\_{\text{dd}} n\_d) - \mu d n\_d + \mu \epsilon n\_{\text{c}}. \end{aligned} \tag{1}$$

The first term on the right hand side of each equation represents the dispersal of the phenotype through diffusion, where *De* and *Dd* are the dispersal rates of each morph. The second term describes the growth rate of the phenotype using a logistic term, this is similar to what is used in Fisher's model [16]. We use *re* and *rd* to represent the growth rate of each morph, *mee* and *mdd* represent the intramorph competition, while *med* and *mde* represent the intermorph competition. The third and fourth terms represent a linear mutation between the phenotypes at mutation rates of *μe* and *μd*, where *μ*, *e* and *d* are constants. Note that we slightly modify the model in the work by the authors of [3] here by replacing the parameters *μe*, *μ<sup>d</sup>* in the work by the authors of [3] with *μe* and *μd* to enable dependence on mutation to be investigated by variation of the single parameter *μ*. It is assumed that all parameters of the system are positive real numbers.

As in the work by the authors of [3], we suppose a basic trade-off between dispersal and growth, namely, that the establisher phenotype has the larger growth rate, while the disperser phenotype has the larger dispersal rate,

$$r\_d > r\_{d\prime} \quad D\_d > D\_{\varepsilon}.\tag{2}$$

While trade-off (2) is not needed either in the proof of linear determinacy or in some of our results on the dependence of spreading speed on mutation rate, we will make use of (2) to discuss parameter-dependent options for the vanishing-mutation limit of spreading speeds in Section 3, and in Section 4, where we characterize the composition of the leading edge of solutions of (1). Further discussion on interesting possible implications of dispersal–growth trade-offs for this model is presented in [12]. Following classical competition theory [19], we suppose throughout that the intramorph competition is greater than the intermorph competition,

$$m\_{\rm dd} \; \geqslant \; m\_{\rm cd}, \quad m\_{\rm cc} \; \geqslant \; m\_{\rm dc}. \tag{3}$$

We also have in mind throughout that the mutation rate *μ* is relatively small in comparison to the other parameters, to remain biologically realistic [14].

Kolmogorov, Petrovskii and Piskunov [15] studied the existence of monotonic travelling wave solutions of the scalar form of the equation

$$\frac{\partial \mathfrak{u}}{\partial t} = A \frac{\partial^2 \mathfrak{u}}{\partial \mathbf{x}^2} + f(\mathfrak{u}). \tag{4}$$

Throughout this work, we will consider travelling wave solutions to be solutions of the Equation (4) of the form *<sup>u</sup>*(*x*, *<sup>t</sup>*) = *<sup>w</sup>*(*<sup>x</sup>* <sup>−</sup> *ct*), where *<sup>w</sup>* : <sup>R</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup>*<sup>n</sup>* is called the wave profile and *<sup>c</sup>* <sup>∈</sup> <sup>R</sup> is the speed of the wave. Kolmogorov, Petrovskii and Piskunov studied the case when *n* = 1, *A* = *d* and *f*(*u*) = *ru*(1 − *u*), proposed by Fisher [16], and proved there is a continuum of values of *c* for which a monotonic travelling wave solution exists, specifically if *c* ≥ *c*∗, where *c*∗ = 2 √ *rd* is the minimal travelling wave speed, as well as establishing stability properties of the minimal-speed front. Aronson and Weinberger [17] further studied this system and characterised *c*∗ as a spreading speed. These results were extended to cooperative systems of equations for a suitable class of nonlinearities *f* by Volpert, Volpert and Volpert [20].

The system (1) is of the form (4) if we let *<sup>u</sup>* = (*ne*, *nd*)*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>2, *<sup>A</sup>* be a diagonal matrix containing the dispersal rates,

$$A = \begin{pmatrix} D\_d & 0 \\ 0 & D\_d \end{pmatrix},\tag{5}$$

and *f* be a nonlinear function containing the growth, competition and mutation terms,

$$f(n\_{\varepsilon}, n\_d) = \begin{pmatrix} r\_{\varepsilon} n\_{\varepsilon} (1 - m\_{\varepsilon \varepsilon} n\_{\varepsilon} - m\_{\varepsilon d} n\_d) - \mu n\_{\varepsilon} + \mu d n\_d \\ r\_d n\_d (1 - m\_{d\varepsilon} n\_{\varepsilon} - m\_{dd} n\_d) + \mu n\_{\varepsilon} - \mu d n\_d \end{pmatrix}. \tag{6}$$

In the following, we will use the notation *u* > *v* to denote that the *i*th component of each vector satisfies *ui* <sup>&</sup>gt; *vi*, for each *<sup>i</sup>*, similarly for *<sup>u</sup>* <sup>&</sup>lt; *<sup>v</sup>*, *<sup>u</sup>* <sup>≥</sup> *<sup>v</sup>* and *<sup>u</sup>* <sup>≤</sup> *<sup>v</sup>*. We say that *<sup>u</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* is positive if *u* > 0. The notation *u* ∈ (*a*, *b*] denotes that the *i*th component of each vector satisfies the inequality *ai* < *ui* ≤ *bi* for each *i*.

Elliott and Cornell investigated numerically the effect of varying the parameters on the spreading speed of the system and interestingly found that, for certain values of growth and dispersal rate, the system would spread faster in the presence of both phenotypes than just one phenotype would spread in the absence of mutation [3]. They predict the spreading speed obtained for each set of parameters in the limit of small mutation, using the front propagation method of van Saarloos [21], making the assumption that the spreading speed of system (1) is linearly determinate in order to do so. As *μ* → 0 the three possible limiting speeds are

$$v\_{\varepsilon} = 2\sqrt{r\_{\varepsilon}D\_{\varepsilon}} \qquad v\_{d} = 2\sqrt{r\_{d}D\_{d}} \qquad v\_{f} = \frac{|r\_{\varepsilon}D\_{d} - r\_{d}D\_{\varepsilon}|}{\sqrt{(r\_{\varepsilon} - r\_{d})(D\_{d} - D\_{\varepsilon})}} \tag{7}$$

Condition (2) is enough to ensure that *vf* exists and is faster than *ve* and *vd*, which are the spreading speeds of the two Fisher–KPP equations that would be satisfied by each phenotype in isolation. The faster speed *vf* is predicted for parameters in the region of the positive quadrant of (*rd*/*re*, *De*/*Dd*)-space, which satisfies the inequalities

$$\frac{D\_d}{D\_t} + \frac{r\_d}{r\_t} > 2, \quad \frac{D\_c}{D\_d} + \frac{r\_c}{r\_d} > 2,\tag{8}$$

represented by the shaded area in Figure 1.

**Figure 1.** Parameter regions showing when the faster invasion speed observed by the authors of [3] occurs. In the upper left unshaded region the solution travels at the speed at which the establisher travels without competition. In the lower right unshaded region the solution travels at the speed at which the disperser travels without competition. In the shaded region the faster spreading speed is observed.

The spreading speed of the system is said to be linearly determinate if it is the same as the spreading speed of the system obtained when (1) is linearised about the (0, 0) equilibrium, namely,

$$\begin{split} \frac{\partial n\_{\ell}}{\partial t} &= D\_{\ell} \frac{\partial^{2} n\_{\ell}}{\partial \mathbf{x}^{2}} + (r\_{\ell} - \mu \mathbf{e}) n\_{\ell} + \mu d n\_{d} \\ \frac{\partial n\_{d}}{\partial t} &= D\_{d} \frac{\partial^{2} n\_{d}}{\partial \mathbf{x}^{2}} + (r\_{d} - \mu d) n\_{d} + \mu e n\_{\ell} . \end{split} \tag{9}$$

This assumption is one that is suggested by the numerical studies in [3], but is not always true even in the scalar case [21–24]. Stokes [25] calls the minimal wave speed *c*∗ "pulled" if it is equal to the linearised spreading speed, that is, the speed of the front is determined by the individuals at the leading edge. Similarly the minimal wave speed is said to be "pushed" if its speed is greater than the linearised spreading speed, in this case the speed is determined by individuals behind the leading edge. Typically there are qualitative differences in wave behaviour depending on whether the wave is pushed or pulled, e.g., stability in the scalar case is discussed in the work by the authors of [26].

In the case of systems, most sufficient conditions for linear determinacy require a cooperative assumption on the system. A system is cooperative when the off-diagonal elements of the Jacobian matrix of *f* are always non-negative, i.e.,

$$\frac{\partial f\_i(u)}{\partial u\_j} \ge 0, \quad \text{if } i \ne j. \tag{10}$$

In biological terms this would mean that each phenotype benefits from the presence of others. A cooperative system is useful mainly due to the existence of a comparison principle for such systems [27,28], which is useful in particular in the proof of linear determinacy.

**Theorem 1** (Comparison Principle [13], Theorem 3.1)**.** *Let A be a positive definite diagonal matrix. Assume that f is a vector-valued function in* R*<sup>n</sup> that is continuous and piecewise continuously differentiable in* R*,* *and that the underlying system* (4) *is cooperative. Suppose that <sup>u</sup>*(*x*, *<sup>t</sup>*) *and <sup>v</sup>*(*x*, *<sup>t</sup>*) *are bounded on* <sup>R</sup> <sup>×</sup> [0, <sup>∞</sup>) *and satisfy*

$$\frac{\partial u}{\partial t} - A \frac{\partial^2 u}{\partial x^2} - f(u) \le \frac{\partial v}{\partial t} - A \frac{\partial^2 v}{\partial x^2} - f(v)$$

*If u*(*x*, *<sup>t</sup>*0) <sup>≤</sup> *<sup>v</sup>*(*x*, *<sup>t</sup>*0) *for x* <sup>∈</sup> <sup>R</sup>*, then*

$$u(\mathbf{x}, t) \le v(\mathbf{x}, t), \quad \text{for } \mathbf{x} \in \mathbb{R}, \ t \ge t \mathbf{0}.$$

Linear determinacy was shown to hold for some cooperative systems by Lui [27], and the result was extended to more general cooperative systems by Weinberger, Lewis and Li [28], who assumed, in particular, that for any positive eigenvector *q* of *f* (0),

$$f(aq) \le af'(0)q,\quad \text{for } a > 0. \tag{11}$$

This reduces to a condition imposed by Hadeler and Rothe [23] in the scalar case.

Unfortunately, we see from the Jacobian matrix (12) that our system (1) is only partially cooperative,

$$J\_f(n\_\varepsilon, n\_d) = \begin{pmatrix} r\_\varepsilon(1 - 2m\_{\varepsilon\varepsilon}n\_\varepsilon - m\_{\varepsilon d}n\_d) - \mu\varepsilon & \mu d - r\_\varepsilon m\_{\varepsilon d}n\_\varepsilon\\ \mu\varepsilon - r\_d m\_{d\varepsilon}n\_d & r\_d(1 - m\_{d\varepsilon}n\_\varepsilon - 2m\_{\text{d\'d}}n\_d) - \mu d \end{pmatrix} . \tag{12}$$

In fact it is typically only cooperative at small population densities due to the relative smallness of the mutation term (see Figure 2).

**Figure 2.** Solution of the system (1) with parameter values: *De* = 0.3, *Dd* = 1.5, *re* = 1.1, *rd* = 0.2, *mee* = 1.0/1.2, *mdd* = 1.0, *med* = 0.8, *mde* = 0.7, *μe* = 0.001, *μd* = 0.00025. A Heaviside step function was used as initial condition for each component.

However, Girardin [11] has recently established a linear determinacy result for a class of non-cooperative reaction–diffusion systems that includes the case considered here; see Theorem 1.7, together with Theorems 1.5 and 1.6, of Girardin [11], which apply here because when (*ne*, *nd*)=(0, 0), the Jacobian (12) always has at least one positive eigenvalue, which can be seen from arguments similar to those used in the proof of Proposition 3.1 below. An alternative proof of linear determinacy for the particular system (1) is given by Morris [12], Theorems 4.5 and 2.15, using the linear determinacy framework outlined by Wang [13]. This latter result is proved under the assumption of sufficiently small mutation and intermorph competition, and uses the fact that (1) is cooperative at low population densities to trap the nonlinearity *f* between two cooperative nonlinearities, *f* <sup>−</sup> and *f* <sup>+</sup>.

Exploiting this linear determinacy, we answer ecologically important questions pertaining to our system in the case of small mutation rate. We first take advantage of a Perron–Frobenius structure to investigate the effect of mutation on spreading speed, and show that an increase in mutation between morphs results in a decrease in the value of the spreading speed; see Theorem 3. This slowing of the speed of propagation as the mutation rate increases is mathematically related, in fact, to the so-called 'reduction phenomenon', that greater mixing lowers growth, discussed in Altenberg [29]. Secondly we investigate the composition of the leading edge of invasion in the limit of small mutation rate, and demonstrate the effects of dispersal, growth rate and mutation on this composition. As a by-product, we also characterise the vanishing-mutation limit of the spreading speed for three different regimes of diffusion and growth parameters in Theorem 2, which yields, in particular, a rigorous explanation of the parameter-dependent selection criteria for the three possible limiting speeds (7) that were discussed in the work by the authors of [3].

We draw the reader's attention to two further interesting references that tackle questions for systems related to (1). Griette and Raoul [30], motivated by an epidemiological model, studied the existence and properties of travelling waves for a special case of system (1). They assume, in particular, that *De* = *Dd*, of which advantage can be taken to prove an explicit formula for the spreading speed and to characterise the shape of travelling wave solutions, including proving non-monotone behaviour in one phenotype and asymptotic behaviour at ±∞. Clearly this assumption of equal dispersal rates, though realistic for the modelling of wild type and mutant types of a virus in the work by the authors of [30] and extremely useful mathematically, is not reasonable for the dispersal polymorphism that is our focus here. Cantrell, Cosner and Yu [31] study (1) from the perspective not of propagation phenomena but of dynamics on a bounded domain. They provide a detailed study of equilibria, the phase plane and dynamics for a range of different parameters, including various regimes for the competition parameters *mee*, *med*, *med*, *mde*.

The rest of the paper is organised as follows. Section 2 presents preliminary material on equilibria of the system (1) and their relationship to equilibria for the related competition–diffusion system when *μ* = 0. The effect of mutation on spreading speeds is discussed in Section 3, using the characterisation of the spreading speed as the linearly-determined minimal speed of a family of travelling waves, which can be expressed and analysed in terms of Perron–Frobenius matrix theory. Section 4 focusses on the parameter regime, in which the dispersal and growth of both morphs play a role in the vanishing-mutation limit of the spreading speed and derives an expression for the ratio of the morphs in the leading edge of the invasion in this case. Some conclusions and remarks are given in Section 5.

#### **2. Equilibria of the System**

We begin with a brief discussion of the equilibria of the system (1) under our assumption (3) on the competition parameters; see also the work by the authors of [31] for further investigation of equilibria of (1). A much studied competition–diffusion system, similar to (1) but where there is no mutation between phenotypes and both intramorph competition values equal one, is the Lotka–Volterra system of equations [32–35],

$$\begin{split} \frac{\partial n\_{\varepsilon}}{\partial t} &= D\_{\varepsilon} \frac{\partial^{2} n\_{\varepsilon}}{\partial \mathbf{x}^{2}} + r\_{\varepsilon} n\_{\varepsilon} (1 - n\_{\varepsilon} - m\_{\varepsilon d} n\_{d}), \\ \frac{\partial n\_{d}}{\partial t} &= D\_{d} \frac{\partial^{2} n\_{d}}{\partial \mathbf{x}^{2}} + r\_{d} n\_{d} (1 - n\_{d} - m\_{d\varepsilon} n\_{\varepsilon}). \end{split} \tag{13}$$

Lewis, Li and Weinberger [5] note that a coexistence equilibrium for this system (13) exists if, and only if, (1 − *med*)(1 − *mde*) > 0; that is, either when both *med* < 1 and *mde* < 1, or both *med* > 1 and *mde* > 1. Note that the case where *med* < 1 and *mde* < 1 corresponds to the condition (3) in our system (1). A stability analysis shows that this coexistence equilibrium is stable when *med* < 1 and *mde* < 1, and unstable when *med* > 1 and *mde* > 1, where here stability is understood in the sense of stability of the ODE system given by (13) with *De* = *Dd* = 0. It should also be noted that the works by the authors of [5,34] also study the case where one species invades the territory of another, while we, and also the authors of [3,30], are concerned with two morphs of a species invading a previously unoccupied territory.

For our system (1) we can easily see that there only exist two constant equilibria by plotting the nullclines of (6),

$$
\mu r\_{\varepsilon} n\_{\varepsilon} (1 - m\_{\alpha \varepsilon} n\_{\varepsilon} - m\_{\text{cd}} n\_{\text{d}}) - \mu e n\_{\varepsilon} + \mu d n\_{\text{d}} = 0,\tag{14}
$$

$$
\mu\_d n\_d (1 - m\_{dt} n\_\varepsilon - m\_{dd} n\_d) - \mu dn\_d + \mu en\_\varepsilon = 0,\tag{15}
$$

and observe where they intersect. The nullclines confirm for a specific choice of parameters that we only have two non-negative equilibria and that they are the extinction equilibrium and a single coexistence equilibrium (Figure 3a). The nullclines appear as they do in Figure 3a if the parameters satisfy the conditions

$$\frac{\mu d}{r\_c m\_{cd}} < \frac{r\_c - \mu e}{r\_c m\_{ce}}, \qquad \frac{\mu e}{r\_d m\_{de}} < \frac{r\_d - \mu d}{r\_d m\_{dd}},\tag{16}$$

which aligns with our assumption that the mutation is relatively small. We note that it is clearly also possible to deal with cases in which the mutation does not satisfy assumptions (16), however we are only interested in the case of small mutation here.

**Figure 3.** (**a**) Nullclines of Equation (6) and (**b**) Nullclines of Equation (17). Parameter values *De* = 0.3, *Dd* = 1.5, *re* = 1.1, *rd* = 0.2, *mee* = 1.0/1.2, *mdd* = 1.0, *med* = 0.8, *mde* = 0.7, *μe* = 0.01, *μd* = 0.025. Each point at which the nullclines intersect represents an equilibrium. We can see that for this choice of parameters Equation (6) has only two non-negative equilibria, while (17) has four.

However, simply plotting the nullclines of our system does not tell us the stability of each equilibrium, we therefore consider a modified version of (6) without the mutation terms, which we call *g*,

$$\mathbf{g}(n\_{\ell\prime}, n\_d) = \begin{pmatrix} r\_{\ell} n\_{\ell} (1 - m\_{\ell\ell} n\_{\ell} - m\_{\ell d} n\_d) \\ r\_{d} n\_{\ell} (1 - m\_{d\ell} n\_{\ell} - m\_{\ell d} n\_d) \end{pmatrix} . \tag{17}$$

Due to the relative smallness of the mutation terms, we can then introduce them as a perturbation before using the implicit function theorem.

First we evaluate the equilibria of *g*. We can easily see that there are four equilibria by plotting the nullclines,

$$(r\_{\mathfrak{e}} n\_{\mathfrak{e}} (1 - m\_{\mathfrak{e}\mathfrak{e}} n\_{\mathfrak{e}} - m\_{\mathfrak{e}\mathfrak{e}} n\_{\mathfrak{d}}) = 0,\tag{18}$$

$$
\sigma\_d n\_d (1 - m\_{de} n\_d - m\_{dd} n\_d) = 0,\tag{19}
$$

which we do in Figure 3b for certain parameters. The equilibria of (17) consist of an extinction equilibrium (0, 0), two equilibria on the axes where one phenotype is present while the other is extinct, (1/*mee*, 0), (0, 1/*mdd*), and a coexistence equilibrium

$$\left(\frac{m\_{dd} - m\_{cd}}{m\_{cc}m\_{dd} - m\_{cd}m\_{dc}}, \frac{m\_{cc} - m\_{dc}}{m\_{cc}m\_{dd} - m\_{cd}m\_{dc}}\right)$$

which we refer to as (*n*∗ *<sup>e</sup>* , *n*<sup>∗</sup> *<sup>d</sup>*) for simplicity. Note that (*n*<sup>∗</sup> *<sup>e</sup>* , *n*<sup>∗</sup> *<sup>d</sup>*) is a coexistence equilibria due to the condition (3) specified earlier.

The Jacobian of (17) is

$$J\_{\mathcal{S}}(n\_{\varepsilon}, n\_{d}) = \begin{pmatrix} r\_{\varepsilon}(1 - 2m\_{\rm cl}n\_{\varepsilon} - m\_{\rm cl}n\_{d}) & -r\_{\varepsilon}m\_{\rm cl}n\_{\varepsilon} \\ -r\_{d}m\_{\rm cl}n\_{d} & r\_{d}(1 - m\_{\rm cl}n\_{\varepsilon} - 2m\_{\rm id}n\_{d}) \end{pmatrix}. \tag{20}$$

Substituting in values of *ne* and *nd* at each of the equilibria to the trace and determinant of (20) we see that the equilibrium (*n*∗ *<sup>e</sup>* , *n*<sup>∗</sup> *<sup>d</sup>*) is stable, while the other three are unstable. Note also that the determinant of (20) is non-zero when evaluated at each of the equilibria of (17).

We now use the implicit function theorem [36] to determine how each equilibrium moves when mutation is introduced to the system (17) as a perturbation. To do so, we suppose that there exists *μ* > 0 such that

$$f(n\_{\mathfrak{e}}, n\_d) = g(n\_{\mathfrak{e}}, n\_d) + \mu M \binom{n\_{\mathfrak{e}}}{n\_d} \tag{21}$$

where *g* is defined in (17) above, *μ* is a non-negative scalar parameter which we use to vary the mutation and *M* is the matrix of mutation coefficients

$$M = \begin{pmatrix} -\varepsilon & d \\ \varepsilon & -d \end{pmatrix}. \tag{22}$$

The equilibria for our original system (1) satisfy *f*(*ne*, *nd*) = 0, where *f* is the nonlinearity (6), so that

$$
\mathcal{g}(n\_{\varepsilon}, n\_d) + \mu \mathcal{M} \left( \begin{array}{c} n\_{\varepsilon} \\ n\_d \end{array} \right) = 0. \tag{23}
$$

As a consequence of the implicit function theorem, in a neighbourhood of *μ* = 0 and an equilibrium (*n*¯*e*, *n*¯ *<sup>d</sup>*) of *g*, there is a unique solution of (23) which is a continuously differentiable function of *μ*, say *h*(*n*¯*e*,*n*¯ *<sup>d</sup>*)(*μ*), provided *g* is invertible at (*n*¯*e*, *n*¯ *<sup>d</sup>*). This ensures we can differentiate (23) in order to obtain an expression describing how an equilibrium (*n*¯*e*, *n*¯ *<sup>d</sup>*)*<sup>T</sup>* is perturbed upon the introduction of mutation *μ*. Since the determinant of the Jacobian matrix *Jg* is not equal to zero at any of the equilibria, we may invert *Jg* and obtain the expression

$$\Theta(\vec{n}\_{\epsilon}, \vec{n}\_{d}) := \frac{\mathbf{d}}{\mathbf{d}\mu} h\_{(\mathfrak{n}\_{\epsilon}, \mathfrak{n}\_{d})}(\mu) \bigg|\_{\mu=0} = -J\_{\mathcal{S}}(\vec{n}\_{\epsilon}, \vec{n}\_{d})^{-1} M \left( \begin{array}{c} \vec{n}\_{\epsilon} \\ \vec{n}\_{d} \end{array} \right). \tag{24}$$

Clearly the extinction equilibrium (0, 0) remains at (0, 0), and the implicit function theorem ensures the local uniqueness of this equilibrium for small *μ* > 0. Evaluating (24) at each of the other equilibria of *g*, we see that the equilibrium (1/*mee*, 0) is perturbed into the lower right quadrant, because

$$\Theta\left(\frac{1}{m\_{\rm ce}},0\right) = \frac{\mu e}{r\_{\rm c}r\_{d}(m\_{\rm ce} - m\_{\rm de})} \left(\frac{r\_{\rm c}m\_{\rm cd}}{m\_{\rm ce}} - \frac{r\_{\rm d}}{m\_{\rm ce}}(m\_{\rm ce} - m\_{\rm de})\_{\rm \prime} - r\_{\rm c}\right)^{T}.\tag{25}$$

Note that the term (*mee* − *mde*) is positive due to the condition (3). Similarly, the equilibrium (0, 1/*mdd*) is perturbed into the upper left quadrant, since

$$\Theta\left(0,\frac{1}{m\_{dd}}\right) = \frac{\mu d}{r\_{\varepsilon}r\_{d}(m\_{dd} - m\_{ed})} \left(-r\_{d}, \frac{r\_{d}m\_{de}}{m\_{dd}} - \frac{r\_{\varepsilon}}{m\_{dd}}(m\_{dd} - m\_{ed})\right)^{T}.\tag{26}$$

Finally, the coexistence equilibrium is perturbed a small amount in a direction which is dependant on the parameters of the system,

$$\Theta\left(n\_{\varepsilon}^{\ast}, n\_{d}^{\ast}\right) = \begin{pmatrix} r\_{\varepsilon} n\_{\varepsilon}^{\ast} n\_{d}^{\ast} \left[ \mu e (m\_{dd} - m\_{\rm ed}) - \mu d (m\_{\rm d\varepsilon} - m\_{\rm dc}) \right] \\ r\_{d} n\_{\varepsilon}^{\ast} n\_{d}^{\ast} \left[ \mu d (m\_{\rm ce} - m\_{\rm dd}) - \mu e (m\_{\rm dd} - m\_{\rm ed}) \right] \end{pmatrix}. \tag{27}$$

Moreover, since the Jacobian is a continuous function of *μ*, we know that for small *μ* = 0, the stability of each of the equilibria remains the same as when *μ* = 0. Therefore, by introducing a small amount of mutation to our system, we are left with two non-negative equilibria: an unstable extinction state (0, 0) and a stable coexistence state (*n*∗ *<sup>e</sup>* , *n*<sup>∗</sup> *d*).

#### **3. The Role of Mutation in Spreading Speeds**

In this and the following section, we derive predictions about the spreading of species modelled by (1) by exploiting the linear determinacy of the system together with the fact that the spreading speed can be characterised using travelling waves.

We being by deriving a *μ*-dependent expression for the minimal travelling wave speed of the linearisation of (1) about the origin. If the general reaction–diffusion system (4) admits a travelling wave solution *u*(*x*, *t*) = *w*(*x* − *ct*), then by substituting *w*(*x* − *ct*) in to the general form (4) we may write the system in the form of a travelling wave equation:

$$A w''(\xi) + c w'(\xi) + g(w(\xi)) + \mu M w(\xi) = 0. \tag{28}$$

We now have an ordinary differential equation in the single variable *ξ* = *x* − *ct*, the linearisation of which about the origin is

$$Aw''(\xi) + cw'(\xi) + g'(0)w(\xi) + \mu Mw(\xi) = 0. \tag{29}$$

Further, if we substitute the ansatz solution *w*(*ξ*) = *e*−*βξq* into (29), we obtain the eigenvalue problem

$$\left(\beta A + \beta^{-1}(\lg'(0) + \mu M)\right)q = cq,\tag{30}$$

where *β* > 0 is the spatial decay and *q* > 0 denotes the phenotypic distribution at the leading edge. We define the matrix on the left hand side to be

$$H\_{\beta,\mu} := \beta A + \beta^{-1}(\mathcal{g}'(0) + \mu M) = \beta A + \beta^{-1}f'(0). \tag{31}$$

For *μ* > 0 the matrix *Hβ*,*<sup>μ</sup>* has strictly positive off-diagonal elements and therefore by the Perron–Frobenius theorem has a Perron–Frobenius eigenvalue, which is plotted in Figure 4 as a function of *β*. This Perron–Frobenius eigenvalue, which we denote *η*PF *Hβ*,*<sup>μ</sup>* , is the larger of the two real eigenvalues of *Hβ*,*<sup>μ</sup>* and has a one-dimensional eigenspace spanned by a positive eigenvector. Further, this Perron–Frobenius eigenvalue is positive for every *β*, *μ* > 0.

**Figure 4.** Perron–Frobenius eigenvalue of *Hβ*, which exists when *μ* > 0. Parameters used are *De* = 0.3, *Dd* = 1.5, *re* = 1.1, *rd* = 0.2, *e* = 0.001, *d* = 0.00025, *μ* = 1.

**Proposition 1.** *The Perron–Frobenius eigenvalue of the matrix βA* + *β*−<sup>1</sup> *f* (0) *is strictly positive for all μ* > 0*, β* > 0

**Proof.** Multiplying the matrix *βA* + *β*−<sup>1</sup> *f* (0) on the right by (*d*,*e*)*<sup>T</sup>* we obtain

$$(\beta A + \beta^{-1}f'(0))\left(\begin{array}{c} d\\e \end{array}\right) = \left(\begin{array}{c} \beta D\_{\mathcal{C}}d + \beta^{-1}r\_{\mathcal{C}}d\\\beta D\_{\mathcal{C}}e + \beta^{-1}r\_{\mathcal{C}}e \end{array}\right) > 0. \tag{32}$$

By Corollary 1.6 of Crooks [37], this implies

$$
\eta\_{\rm PF}(\beta A + \beta^{-1} f'(0)) > 0. \tag{33}
$$

**Definition 1.** *The minimal speed of the travelling wave solution w*(*ξ*) = *e*−*βξq of* (28) *for a given μ* > 0 *is given by*

$$\eta\left(\mu\right) := \inf\_{\beta > 0} \eta\_{\text{PF}}\left(H\_{\beta,\mu}\right). \tag{34}$$

*We define β*(*μ*) *to be the value of β at which η*(*μ*) *is attained.*

Note that it follows from Gershgorin's Circle Theorem that *η*PF(*Hβ*,*μ*) → ∞ as *β* → 0 and *β* → ∞, which, together with Lemma 1.1 (5)–(6) of Wang [13] and the fact that *β* → *<sup>η</sup>PF*(*β*2*<sup>A</sup>* + *<sup>f</sup>* (0)) is a strictly convex function of *β* by Lemma 3.7 of Crooks [37], implies that there exists a unique value *β* = *β*(*μ*) at which inf*β*><sup>0</sup> *η*PF(*Hβ*,*μ*) is attained.

In the case *μ* = 0, the matrix *Hβ*,0 is diagonal, and therefore does not have a Perron–Frobenius eigenvalue. Instead the matrix *Hβ*,0 has the two eigenvalues

$$
\beta D\_d + \beta^{-1} r\_{d\prime} \quad \text{and} \quad \beta D\_e + \beta^{-1} r\_{e\prime} \tag{35}
$$

which we plot in Figure 5a as functions of *β* (for certain parameters). The minimum of the first curve in (35) is 2√*rdDd*, which we note is the Fisher speed of the disperser *vd*, and occurs at the *<sup>β</sup>*-value

$$
\beta\_d = \sqrt{r\_d / D\_d}.\tag{36}
$$

Similarly, the minimum of the second curve in (35) is the Fisher speed of the disperser *ve* = 2 <sup>√</sup>*reDe*, and occurs at the *β*-value

$$
\beta\_{\varepsilon} = \sqrt{r\_{\varepsilon}/D\_{\varepsilon}}.\tag{37}
$$

The third point of interest is the point at which the the two curves in (35) cross, which is *vf* , and occurs at the *β*-value

$$
\beta\_f = \frac{\sqrt{r\_c - r\_d}}{\sqrt{D\_d - D\_c}}.\tag{38}
$$

By analogy with the Perron–Frobenius eigenvalue we consider the maximum of these two eigenvalues for each *β*, which we plot in Figure 5b. Note here the similarity to the plot of the Perron–Frobenius eigenvalue of *Hβ*,*μ*, seen in Figure 4.

**Figure 5.** (**a**) Eigenvalues of *H<sup>β</sup>* in the case of *μ* = 0 when *De* = 0.3, *Dd* = 1.5, *re* = 1.1, *rd* = 0.2. (**b**) *η* corresponding to each *β* for *μ* = 0.

**Definition 2.** *We define the minimum over β of the maximum function of eigenvalues defined in* (35) *to be η*0*. Likewise, we denote the value of β at which this minimum is obtained by β*∗*.*

Note that, in Figure 5, the minimum over *β* of the maximum of the eigenvalues (35) is the point at which they both meet, and therefore in this case *η*<sup>0</sup> = *vf* and *β*<sup>∗</sup> = *β <sup>f</sup>* .

However there are also regions of parameters in which the eigenvalues (35) meet in such a way that this is not the case. The condition that the minima of the two eigenvalues lie on either side of the point at which they meet, which is equivalent to the minimum value of *β* of the maximum of the eigenvalues being at the crossing point, as in Figure 5, is

$$
\sqrt{\frac{r\_d}{D\_d}} < \sqrt{\frac{r\_c - r\_d}{D\_d - D\_c}} < \sqrt{\frac{r\_c}{D\_c}}\tag{39}
$$

which implies the condition

$$\frac{r\_c}{r\_d} + \frac{D\_c}{D\_d} > 2 \quad \text{and} \quad \frac{r\_d}{r\_c} + \frac{D\_d}{D\_c} > 2,\tag{40}$$

imposed by Elliott and Cornell [3] ensuring that the faster spreading speed *vf* (7) is obtained in the limit as *μ* → 0. Further, when condition (40) is satisfied, *Hβ*∗,0 is a multiple of the identity and has a repeated eigenvalue with a two-dimensional eigenspace. We consider this case, where *vf* is obtained as the limiting speed as *μ* → 0, to be of most biological interest, as unlike *ve* and *vd*, *vf* is then dependent on the traits of both morphs and occurs as a result of polymorphism.

In the case that the condition (40) is not satisfied we would expect one of the individual speeds be selected. For example if the dispersal and growth rates of each phenotype satisfy

$$\frac{r\_\varepsilon}{r\_d} + \frac{D\_\varepsilon}{D\_d} < 2 \quad \text{and} \quad \frac{r\_d}{r\_\varepsilon} + \frac{D\_d}{D\_\varepsilon} > 2,\tag{41}$$

then the individual speed of phenotype *d* is selected. In this case the minimum over *β* of the maximum of the eigenvalues (35) is no longer the crossing point of the eigenvalues, and is instead the minimum of the eigenvalue curve *βDd* + *β*−<sup>1</sup>*rd*, therefore *η*<sup>0</sup> = *vd* and *β*<sup>∗</sup> = *βd*. Similarly, in the case

$$\frac{r\_\varepsilon}{r\_d} + \frac{D\_\varepsilon}{D\_d} > 2 \quad \text{and} \quad \frac{r\_d}{r\_\varepsilon} + \frac{D\_d}{D\_\varepsilon} < 2,\tag{42}$$

where the individual speed of phenotype *e* is selected, we see that *η*<sup>0</sup> = *ve* and *β*<sup>∗</sup> = *βe*. It is not possible for both inequalities in (40) to be reversed while the condition (2) holds.

**Proposition 2.** *If* (40) *holds, then η*<sup>0</sup> = *vf and β*<sup>∗</sup> = *β <sup>f</sup> . If* (41) *holds, then η*<sup>0</sup> = *vd and β*<sup>∗</sup> = *βd. If* (42) *holds, then η*<sup>0</sup> = *ve and β*<sup>∗</sup> = *βe.*

The values *η*<sup>0</sup> and *β*∗, defined in Definition 2 using the eigenvalues of the matrices *Hβ*,0, are shown in Morris [12] Theorems 5.6–5.11 to be the limits of *η*(*μ*) and *β*(*μ*) as the mutation rate *μ* → 0. In light of the linear determinancy of the spreading speed established in the work by the authors of [11,12], this yields a rigorous characterisation of the limit of the spreading speed as *μ* → 0, which shows, in particular, that the limiting speeds (7) predicted by Elliott and Cornell [3] using the front propagation method of van Saarloos [21] are indeed what is obtained in the limit of small mutation. Here we summarise the results and refer to Theorems 5.6–5.11 of Morris [12] for details of the proofs.

**Theorem 2.** *In all three of the parameter regions* (40)*–*(42)*,*

$$\lim\_{\mu \to 0} \eta(\mu) = \eta\_{0\prime} \quad \lim\_{\mu \to 0} \beta(\mu) = \beta^\* \mu$$

*where η*0, *β*∗ *are as defined in Definition 2, that is,*

*(i) if* (40) *holds, then*

$$\lim\_{\mu \to 0} \eta(\mu) = \frac{|r\_{\mathfrak{e}} D\_d - r\_d D\_{\mathfrak{e}}|}{\sqrt{(r\_{\mathfrak{e}} - r\_d)(D\_d - D\_{\mathfrak{e}})}}, \qquad \lim\_{\mu \to 0} \beta(\mu) = \frac{\sqrt{r\_{\mathfrak{e}} - r\_d}}{\sqrt{D\_d - D\_{\mathfrak{e}}}}.$$

*(ii) if* (41) *holds, then*

$$\lim\_{\mu \to 0} \eta(\mu) = 2\sqrt{r\_d D\_d} \qquad \lim\_{\mu \to 0} \beta(\mu) = \sqrt{r\_d/D\_d}.$$

*(iii) if* (42) *holds, then*

$$\lim\_{\mu \to 0} \eta(\mu) = 2\sqrt{r\_{\varepsilon}D\_{\varepsilon}} \qquad \lim\_{\mu \to 0} \beta(\mu) = \sqrt{r\_{\varepsilon}/D\_{\varepsilon}}.$$

Note that Tang and Fife [38] established the existence of travelling fronts for (1) for all speeds greater than or equal to max{2 <sup>√</sup>*reDe*, 2√*rdDd*}. Therefore, interestingly, in the case when (2) and (40) are satisfied, the limit as *μ* → 0 of the minimal front speed *η*(*μ*) is strictly larger than the minimal front speed for (1) in the absence of mutation.

Since calculating the spreading speed *η*(*μ*) involves minimizing in *β* and the diffusion matrix *A* is not a multiple of the identity, finding an explicit expression for *η*(*μ*) is not very tractable. However, by adapting an argument from the work by the author of [37], we show that *η*(*μ*) is a nonincreasing function of the mutation rate *μ*. A key ingredient of the proof is the following classical result of Cohen [39] on convexity properties of Perron–Frobenius eigenvalues.

**Lemma 1** (Cohen [39])**.** *Let <sup>P</sup>*, *<sup>Q</sup>* <sup>∈</sup> <sup>R</sup>*n*×*<sup>n</sup> such that <sup>P</sup> is diagonal and <sup>Q</sup> has positive off-diagonal elements. Then the Perron–Frobenius eigenvalue of P* + *Q is a convex function of P; that is, given diagonal matrices P*<sup>1</sup> *and P*<sup>2</sup> *and* 0 < *α* < 1*,*

$$a\eta\_{\rm PF}(aP\_1 + (1-a)P\_2 + Q) \le a\eta\_{\rm PF}(P\_1 + Q) + (1-a)\eta\_{\rm PF}(P\_2 + Q). \tag{43}$$

This convexity, together with the fact that *ηPF*(*M*) = 0, can now be exploited to obtain the following.

**Theorem 3.** *η*(*μ*) *is a nonincreasing function of μ.*

**Proof.** Take *μ*<sup>0</sup> > 0, define *Z* to be the 2 × 2 zero matrix, and define

$$P := \beta(\mu\_0)^2 A - \beta(\mu\_0)\eta(\mu\_0)I + \lg'(0) \tag{44}$$

Note that *β*(*μ*0) > 0 and

$$
\eta\_{\rm PF}(\beta^2 A - \beta \eta(\mu\_0)I + \mathbf{g}'(0) + \mu\_0 \mathbf{M}) \ge 0, \quad \forall \beta > 0,\tag{45}
$$

$$
\eta\_{\rm PF} (\beta(\mu\_0)^2 - \beta(\mu\_0)\eta(\mu\_0)I + \mathbf{g}'(0) + \mu\_0 M) = 0. \tag{46}
$$

Now for any *μ* > 0, we know that

$$
\eta\_{\rm PF} \left( P + \mu M \right) = \mu \,\eta\_{\rm PF} \left( \frac{1}{\mu} P + M \right),
\tag{47}
$$

and in particular, *<sup>η</sup>*PF(*<sup>P</sup>* <sup>+</sup> *<sup>μ</sup>M*) and *<sup>η</sup>*PF <sup>1</sup> *<sup>μ</sup> P* + *M* have the same sign. Moreover,

$$
\eta\_{\rm PF} \left( \frac{1}{\mu\_0} P + M \right) = 0. \tag{48}
$$

For any *μ* > *μ*0, we can write,

$$\frac{1}{\mu}P = \frac{\mu\_0}{\mu} \left(\frac{1}{\mu\_0}P\right) + \left(1 - \frac{\mu\_0}{\mu}\right)Z.\tag{49}$$

Now in Lemma 1, let *P*/*μ*<sup>0</sup> and *Z* be *P*<sup>1</sup> and *P*2, respectively, and let *α* := *μ*0/*μ*. Then by the convex dependence of *η*PF,

$$\begin{split} \eta\_{\text{PF}}\left(\frac{1}{\mu}P + M\right) &\leq \frac{\mu\_{0}}{\mu} \eta\_{\text{PF}}\left(\frac{1}{\mu\_{0}}P + M\right) + \left(1 - \frac{\mu\_{0}}{\mu}\right) \eta\_{\text{PF}}\left(Z + M\right) \\ &= \frac{\mu\_{0}}{\mu} \eta\_{\text{PF}}\left(\frac{1}{\mu\_{0}}P + M\right) + \left(1 - \frac{\mu\_{0}}{\mu}\right) \eta\_{\text{PF}}\left(M\right) . \end{split} \tag{50}$$

We know from (48) that the first term is equal to zero, and *η*PF(*M*) = 0 implies the second term is also zero. The inequality (50) therefore becomes

$$
\eta\_{\rm PF} \left( \frac{1}{\mu} P + M \right) \le 0, \quad \mu > \mu\_0. \tag{51}
$$

Using (47), this implies

$$
\eta\_{\rm PF} \left( P + \mu M \right) \le 0, \qquad \mu > \mu\_0, \tag{52}
$$

and substituting in the full expression for *P* and dividing by *β*(*μ*0) then yields

$$
\eta\_{\rm PF}(\beta(\mu\_0)A + \beta(\mu\_0)^{-1}(\mathbf{g}'(0) + \mu M)) \le \eta(\mu\_0), \quad \mu > \mu\_0. \tag{53}
$$

Hence

$$\inf\_{\beta>0} \eta\_{\text{FF}}(\beta A + \beta^{-1}(\mathbf{g}'(0) + \mu M)) \le \eta(\mu\_0). \tag{54}$$

Since the expression on the left hand side of (54) is the definition of *η*(*μ*), we have

$$
\eta(\mu) \le \eta(\mu\_0). \tag{55}
$$

Since this holds for any *μ* > *μ*0, we have therefore shown that *η*(*μ*) is a nonincreasing function of *μ*.

In light of the linear determinancy of the spreading speed [11,12], Theorem 3 establishes that the spreading speed for the nonlinear system (1) is a nonincreasing function of the mutation rate *μ*. Interestingly, this is related mathematically to the so-called "reduction phenomenon" discussed by Altenberg [29], which roughly says that, under certain conditions, greater mixing results in lowered growth. On the other hand, our result can be interpreted as showing that greater mixing results in a slower speed of propagation. It is possible, in fact, to make use of Theorem 6(iii) [29] to give a slightly different proof of Proposition 3.

#### **4. Behaviour at the Leading Edge in the Limit** *μ* → 0

Here we study the Perron–Frobenius eigenvector *qβ*(*μ*) corresponding to the eigenvalue *η*(*μ*) in the limit as the mutation rate *μ* → 0, with the aim of determining the ratio of the phenotypes in the leading edge of an invasion. We do this under conditions (2) and (8) on the dispersal and growth parameters, which, by Theorem (2)(i), ensures that the faster speed *vf* is obtained as *μ* → 0 and, as we will see, results in both phenotypes being present in the leading edge. Throughout this section we therefore assume that *re*

$$\frac{r\_\varepsilon}{r\_d} + \frac{D\_\varepsilon}{D\_d} > 2 \quad \text{and} \quad \frac{r\_d}{r\_\varepsilon} + \frac{D\_d}{D\_\varepsilon} > 2,\tag{56}$$

as well as

$$r\_{\mathfrak{e}} > r\_{\mathfrak{d}\prime} \; D\_{\mathfrak{d}} > D\_{\mathfrak{e}\prime} \tag{57}$$

and for convenience in the following, define the quantities

$$a := \left\| \boldsymbol{\upmu}^2 \boldsymbol{D}\_d - \boldsymbol{r}\_d, \quad \boldsymbol{b} := -\boldsymbol{\upbeta}^{\*^2} \boldsymbol{D}\_d + \boldsymbol{r}\_d. \tag{58}$$

Note that a straightforward calculation shows that the condition

$$a > 0 \quad \text{and} \quad b > 0 \tag{59}$$

is equivalent to the condition (56). We consider this parameter regime to be of greatest biological interest, since unlike *ve* and *vd*, the value *vf* is dependent on the traits of both morphs and occurs as a result of polymorphism. Note that it is shown by (Girardin [40], Theorem 1.1) that the eigenvector *qβ*(*μ*), which arises in the explicit solution of the linearised problem (29), does indeed also characterise the asymptotic behaviour of travelling wave solutions of the nonlinear system (1) close to the extinction state.

In the absence of mutation (*μ* = 0), Figure 5 illustrates the *β*-dependence of eigenvalues of the matrix *Hβ*,0 for a choice of dispersal and growth parameters for which conditions (56) and (57) are satisfied. Under these conditions, the minimal travelling wave speed *η*<sup>0</sup> is obtained at the point at which the two eigenvalues of *Hβ*,0 meet, so that

$$
\eta\_0 = \left\pounds^\* D\_d + \left\pounds^{\*^{-1}} r\_d = \left\pounds^\* D\_\varepsilon + \left\pounds^{\*^{-1}} r\_{\varepsilon\_\prime} \right. \right. \qquad \left. \left. \begin{array}{c} \beta^\* = \sqrt{\frac{r\_\varepsilon - r\_d}{D\_d - D\_\varepsilon}}. \end{array} \right. \tag{60}
$$

In this case, *Hβ*∗,0 is a multiple of the identity, so has a repeated eigenvalue with a two-dimensional eigenspace, making it not obvious a priori what ratio between the two components one would expect in the limit of vanishing mutation.

We therefore investigate the Perron–Frobenius eigenvector *q* of *Hβ*(*μ*),*<sup>μ</sup>* in the limit *μ* → 0, restricting attention to the region of parameters in which (56) and (57) hold. Recall from Theorem 2 that lim*μ*→<sup>0</sup> *η*(*μ*) = *η*<sup>0</sup> and lim*μ*→<sup>0</sup> *β*(*μ*) = *β*∗, and assume that the following limits exist.

$$\eta'(0) = \lim\_{\mu \to 0} \frac{\eta(\mu) - \eta\_0}{\mu}, \quad \beta'(0) = \lim\_{\mu \to 0} \frac{\beta(\mu) - \beta^\*}{\mu}. \tag{61}$$

Note that we provide further numerical justification of these assumptions in Figure 6. Rewriting (30) as a system of scalar equations and taking *β* = *β*(*μ*), we obtain

$$\left(\beta(\mu)D\_{\varepsilon} + \frac{(r\_{\varepsilon} - \mu\varepsilon)}{\beta(\mu)} - \eta(\mu)\right) + \frac{\mu d}{\beta(\mu)}\frac{q\_{2}}{q\_{1}} = 0,\tag{62}$$

$$\left(\beta(\mu)D\_d + \frac{(r\_d - \mu d)}{\beta(\mu)} - \eta(\mu)\right)\frac{q\_2}{q\_1} + \frac{\mu c}{\beta(\mu)} = 0,\tag{63}$$

for *μ* > 0. Then adding and subtracting *η*<sup>0</sup> from (62) and (63), and dividing by *μβ*(*μ*)−<sup>1</sup> yields

$$\beta^\* \frac{(\beta D\_\varepsilon + \beta^{-1} r\_\varepsilon) - (\beta^\* D\_\varepsilon + \beta^{\*^{-1}} r\_\varepsilon)}{\mu} - \beta^\* \frac{\eta - \eta\_0}{\mu} - \varepsilon + d \frac{q\_2}{q\_1} = 0 \tag{64}$$

$$\left(\beta^\* \frac{(\beta D\_d + \beta^{-1} r\_d) - (\beta^\* D\_d + \beta^{\*^{-1}} r\_d)}{\mu} - \beta^\* \frac{\eta - \eta\_0}{\mu} - d\right) \frac{q\_2}{q\_1} + c = 0 \tag{65}$$

Now bearing in mind (61) and the fact that *β*<sup>∗</sup> = 0, it follows from (64) that

$$\frac{q\_{0\_2}}{q\_{0\_1}} := \lim\_{\mu \to 0} \frac{q\_2}{q\_1}$$

exists, so that taking the limit *μ* → 0 in (64), (65), we obtain a system of two equations in the three unknowns *q*02/*q*01 , *η* (0) and *β* (0),

$$\beta^\* \left[ D\_\varepsilon \beta'(0) - \frac{r\_\varepsilon \beta'(0)}{\beta^{\*^2}} - \eta'(0) \right] - \varepsilon + d\left(\frac{q\_{0\_2}}{q\_{0\_1}}\right) = 0,\tag{66}$$

$$
\beta^\* \left[ D\_d \beta'(0) - \frac{r\_d \beta'(0)}{\beta^{\*^2}} - \eta'(0) \right] \left( \frac{q\_{0\_2}}{q\_{0\_1}} \right) - d \left( \frac{q\_{0\_2}}{q\_{0\_1}} \right) + d = 0. \tag{67}
$$

To obtain a third equation, recall that the eigenvalues *λ* of *Hβ*,*<sup>μ</sup>* satisfy det(*Hβ*,*<sup>μ</sup>* − *λI*) = 0, so that

$$(\lambda - \beta D\_{\varepsilon} - \beta^{-1}(r\_{\varepsilon} - \mu \varepsilon))(\lambda - \beta D\_{d} - \beta^{-1}(r\_{d} - \mu d)) - \beta^{-2}\mu^{2}\varepsilon d = 0,\tag{68}$$

and further, at *β* = *β*(*μ*), we have

$$\frac{\partial \lambda}{\partial \beta} = 0,\text{ when }\lambda = \eta(\mu). \tag{69}$$

Differentiating (68) with respect to *β* and setting *β* = *β*(*μ*), *λ* = *η*(*μ*) , we then multiply the result by *β*<sup>3</sup> to obtain

$$\begin{split} \left(\beta^2(D\_\varepsilon + D\_d) - (r\_\varepsilon + r\_d - \mu e - \mu d)\right) \beta \, \eta(\mu) - 2\mu^2 ed &= \\ \left(\beta^2 D\_\varepsilon - (r\_\varepsilon - \mu e)\right) \left(\beta^2 D\_d + (r\_d - \mu d)\right) &+ \left(\beta^2 D\_\varepsilon + (r\_\varepsilon - \mu e)(\beta^2 D\_d - (r\_d - \mu d))\right) .\end{split} \tag{70}$$

Then differentiating with respect to *μ* and letting *μ* → 0 yields a third equation relating *β* (0) and *η* (0), namely

$$2\beta^\*(bD\_d - aD\_\varepsilon)\beta'(0) + \eta\_0(a-b)\beta'(0) + \beta^\*(a-b)\eta'(0) = (bd - ac),\tag{71}$$

where *a* and *b* are as defined in (58). Recall from (59) that *a*, *b* are both positive due to the assumption (56) on the parameters *De*, *Dd*,*re* and *rd*.

We thus have a system of three equation (66), (67) and (71), which we solve. Eliminating *η* (0) and *β* (0) gives the explicit expression for *q*02/*q*<sup>01</sup> :

$$q\_{\text{ratio}} := \frac{q\_{0\_2}}{q\_{0\_1}} = \sqrt{\frac{bc}{ad}} = \sqrt{\frac{(r\_\epsilon - \beta^{\*2} D\_\epsilon)e}{(\beta^{\*2} D\_d - r\_d)d}} > 0. \tag{72}$$

Recall that this is the ratio of the two components of the eigenvector of (30) in the limit *μ* → 0, where *q*<sup>2</sup> represents the disperser and *q*<sup>1</sup> the establisher. This quantity therefore predicts the ratio of the two morphs present in the leading edge of the minimal-speed travelling wave solution to our system [40], which gives insight into the expected long-time behaviour of the solution of (1) with Heaviside initial conditions that models the invasion of the two morphs into an empty region. Note that since *q*ratio > 0, both morphs play a role in the leading edge.

The presence of the two positive quantities *a* and *b* in (72) emphasizes that this expression holds only in the zone where (56), which is equivalent to (59), is satisfied, in which case the faster speed *vf* is obtained in the limit *μ* → 0. If instead either condition (41) or condition (42) holds, in which case the *μ* → 0 limits of *β*(*μ*) and *η*(*μ*) are given by Theorem 3.5 (ii) or (iii) and depend either only on *De*,*re* or only on *Dd*,*rd*, then up to normalisation, the *μ* → 0 limit *q*<sup>0</sup> of the eigenvector *qβ*(*μ*) is either (1, 0)*<sup>T</sup>* or (0, 1)*T*, corresponding to *q*ratio = 0 or ∞. This is because the diagonal matrix *Hβ*∗,0 = diag(*β*∗*De* + *re <sup>β</sup>*<sup>∗</sup> , *<sup>β</sup>*<sup>∗</sup>*Dd* <sup>+</sup> *rd <sup>β</sup>*<sup>∗</sup> ) is then no longer a multiple of the identity and the eigenvector *q*<sup>0</sup> satisfies *Hβ*∗,0 *q*<sup>0</sup> = *η*<sup>0</sup> *q*0. So in these parameter zones, one of the two phenotypes will dominate the leading edge of the front in the limit *μ* → 0, in contrast to the zone when (56) holds, when (72) is valid and both components play a role.

#### *4.1. Discussion of Leading Edge Behaviour*

In order to investigate the effects of parameters on the proportion of each morph in the leading edge, note that using the substitutions *r* = *rd*/*re*, *D* = *De*/*Dd*, *m* = *e*/*d*, we can rewrite (72) as a function of only three variables, the ratio of the growth rates *r*, the ratio of the dispersal rates *D*, and the ratio of the mutation rates *m*, namely,

$$q\_{\text{ratio}} = \sqrt{\left(\frac{2D - rD - 1}{2r - rD - 1}\right)m}.\tag{73}$$

This is both biologically interesting and mathematically elegant.

**Figure 6.** Numerical solutions of (**a**) *η* (*μ*), (**b**) *β* (*μ*) demonstrating convergence as *μ* → 0, obtained using Mathematica. Parameter values are *re* = 1.1, *rd* = 0.2, *De* = 0.3, *Dd* = 1.5, *e* = 0.001 and *d* = 0.00025.

It is instructive and biologically relevant to comment on the effect on *q*ratio of varying parameters *r*, *D*, and *m*. Figure 7 illustrates how the ratio of growth, dispersal and mutation rates affects the proportion of each morph in the leading edge of the invasion wave.

We see in Figure 7a that as *re* increases (or *rd* decreases) there is an increase in the establisher morph in the leading edge. However, as *q*ratio does not fall below 1, there is always a larger amount of the disperser morph in the leading edge. As we increase *re* to infinity (or decrease *rd* to 0) *q*ratio tends to (<sup>1</sup> <sup>−</sup> 2/*D*)*m*. Observe also that in (i), *<sup>q</sup>*ratio blows up at the value of *<sup>r</sup>* <sup>=</sup> *rd*/*re* at which, for the fixed value of *D* = *De*/*Dd* in the plot, our assumption that *a* > 0 and *b* > 0 stops holding. This corresponds to the (normalised) eigenvector *q*<sup>0</sup> approaching the eigenvector (0, 1)*<sup>T</sup>* and the parameters *D*,*r* entering the zone where (41) holds and lim*μ*→<sup>0</sup> *η*(*μ*) is the Fisher speed of the disperser morph *vd* instead of *vf* .

**Figure 7.** Parameter sweeps of Equation (73). We fix *D* and *m* in (**a**), *r* and *m* in (**b**), and *r* and *D* in (**c**). When fixed, parameters take the values *r* = 0.2/1.1, *D* = 0.3/1.5, *m* = 0.001/0.00025.

Figure 7b tells us that as we increase *Dd* (or decrease *De*), there is an increase in the disperser morph in the leading edge. If *De* and *Dd* are close enough in value, it is possible that the establisher morph outnumbers the disperser in the leading edge, which we see in the region in which *q*ratio drops below one. As we increase *Dd* to infinity (or decrease *De* to 0) *<sup>q</sup>*ratio tends to *m*/(<sup>1</sup> <sup>−</sup> <sup>2</sup>*r*). In the case of (ii), the ratio *q*ratio becomes zero at the value of *D* = *De*/*Dd* at which, for the fixed value of *r* = *rd*/*re* in the plot, our assumption that *a* > 0 and *b* > 0 no longer holds. This corresponds to the (normalised) eigenvector *q*<sup>0</sup> approaching the eigenvector (1, 0)*<sup>T</sup>* and the parameters *D*,*r* entering the zone where (42) holds and lim*μ*→<sup>0</sup> *η*(*μ*) is the Fisher speed of the establisher morph *ve* instead of *vf* .

In Figure 7c we see that if the mutation rate of a morph is increased it results in a decrease in the density of that morph in the leading edge. Since the ratio *m* is not restricted to a particular range, like *r* and *D* due to (2), we see that the mutation rate is also able to change which morph is more prevalent in the leading edge. If we fix all parameters but the mutation rates in (73) we obtain *<sup>q</sup>*ratio <sup>=</sup> <sup>√</sup>*km*, where *<sup>k</sup>* <sup>&</sup>gt; 0 is a constant. Thus *<sup>q</sup>*ratio goes to zero as *<sup>m</sup>* grows smaller and blows up as *<sup>m</sup>* tends to infinity. These provide us with interesting predictions which could be tested experimentally with real biological systems.

Note that, of course, the results of the works by the authors of [11,12] on linear determinacy and existence of travelling waves are valid for all parameter choices *De*, *Dd*,*re*,*rd*. The fact that the expression *q*ratio only holds in some parameter zones is just an expression of the fact that only for certain parameters does the eigenvector that describes the leading edge of the front, which is strictly positive whenever *μ* > 0, tend to a strictly positive eigenvector as *μ* → 0.

Recall that we made assumptions (61) on the existence of *η* (0) and *β* (0), for which we now provide some justification by numerically solving Equations (62), (68) and (70) as a system of three equations in *q*2/*q*1, *η*(*μ*) and *β*(*μ*). Differentiating the latter two with respect to *μ* we plot *η* (*μ*), *β* (*μ*) and *q*2/*q*<sup>1</sup> as functions of *μ* in Figure 6 demonstrating the existence of *η* (0), *β* (0) and *q*ratio for this set of parameters. Note that as *μ* increases we see an increase in the proportion of the phenotype *nd*, which is expected as for this set of parameters *e* > *d*.

Finally, we also plot *β*(*μ*) and *η*(*μ*) as functions of *μ* in Figure 8. As *μ* → 0 we see *β*(*μ*) and *η*(*μ*) converge to the values *β*∗ and *η*0, respectively. Further, as we increase *μ* we observe a decrease in the minimal spreading speed *η*(*μ*), this is consistent with our result in Section 4.1.

**Figure 8.** Numerical solutions of (**a**) *β*(*μ*) and (**b**) *η*(*μ*) demonstrating convergence to *β*<sup>∗</sup> and *η*0, respectively, as *μ* → 0, obtained using Mathematica. The behaviour of *η*(*μ*) is consistent with our result that *η*(*μ*) is a nonincreasing function of *μ*. Parameter values: *re* = 1.1, *rd* = 0.2, *De* = 0.3, *Dd* = 1.5, *e* = 0.001 and *d* = 0.00025.

#### **5. Conclusions and Remarks**

This article exploits linear determinacy of the competition–diffusion-mutation system (1) to establish ecologically-motivated results about the invasion properties of a species with two morphs. These results constitute theoretical predictions of the model, of which the most significant are


It would, of course, be valuable to compare these predictions with experimental or empirical data. Note that on the one hand, for the parameter regime (40), so-called 'anomalous spreading' in [3], the spreading speed in the vanishing mutation limit is faster than either phenotype would spread in isolation. On the other hand, once mutation is present, the spreading speed is a nonincreasing function of mutation. Although these predictions are rigorous consequences of the model, their juxtaposition is slightly counter-intuitive from an ecological point of view.

Activity on propagation phenomena in models incorporating mutation has increased markedly in recent years (see the works by the authors of [10,11,30,40] and many others) and is producing a growing body of interesting results and open questions. A natural progression of our work would be a detailed analysis of systems for *N* morphs, in particular, dependence of spreading speeds and the composition of the morphs in the leading edge on parameters. The methods developed in Sections 3 and 4 have the clear potential to yield results in the *N* morph case. Also interesting is the role of functional trade-offs in determining spreading speeds in competition–diffusion mutation models, some first results on which are established in [12]. More wide-ranging questions include how best to model mutation (for instance, should the rate *μ* be replaced by a density-dependent rate), whether continuous or discrete trait models are most appropriate in a given setting, and how to pass rigorously from a discrete setting with a large number of traits to the continuous-trait setting.

**Author Contributions:** All authors contributed equally and significantly in writing this article.

**Funding:** This research received no external funding.

**Acknowledgments:** We would like to acknowledge Natalia Petrovskaya at the University of Birmingham and Sergei Petrovskii at the University of Leicester for their feedback and for organising the META Mathematical Ecology: Theory and Applications series of events. We would also like to thank Léo Girardin at Université Paris-Sud, Stephen Cornell and Vincent Keenan at the University of Liverpool, Jason Matthiopoulos at the University of Glasgow, and Jonathan Potts at the University of Sheffield for their insight and several interesting discussions. The first author would like to acknowledge Swansea University and the College of Science at Swansea University for financial support in the form of a PhD studentship held in the Swansea Centre for Biomathematics. We are also grateful to the two reviewers for useful remarks that have helped to improve the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
