Toxin-mediated competition in weakly motile bacteria

TL;DR

This paper models toxin-mediated competition between bacteria, showing how toxin strength and motility influence colony interactions, inhibition zones, and invasion dynamics through analytical and numerical methods.

Contribution

It extends classic competition models to include toxin effects, deriving analytical solutions for non-motile bacteria and analyzing wave behavior in weakly motile bacteria.

Findings

Inhibition zones form when toxin strength exceeds a threshold.

Wave speed depends on toxin strength and competition, with producer colonies potentially outcompeting susceptibles.

Higher bacterial motility requires stronger toxins for effective inhibition.

Abstract

Many bacterial species produce toxins that inhibit their competitors. We model this phenomenon by extending classic two-species Lotka-Volterra competition in one spatial dimension to incorporate toxin production by one species. Considering solutions comprising two adjacent single-species colonies, we show how the toxin inhibits the susceptible species near the interface between the two colonies. Moreover, a sufficiently effective toxin inhibits the susceptible species to such a degree that an `inhibition zone' is formed separating the two colonies. In the special case of truly non-motile bacteria, i.e. with zero bacterial diffusivity, we derive analytical expressions describing the bacterial distributions and size of the inhibition zone. In the more general case of weakly motile bacteria, i.e. small bacterial diffusivity, these two-colony solutions become travelling waves. We employ numerical methods to show that the wavespeed is dependent upon both interspecific competition and toxin strength; precisely which colony expands at the expense of the other depends upon the choice of parameter values. In particular, a sufficiently effective toxin allows the producer to expand at the expense of the susceptible, with a wavespeed magnitude that is bounded above as the toxin strength increases. This asymptotic wavespeed is independent of interspecific competition and due to the formation of the inhibition zone; when the colonies are thus separated, there is no longer direct competition between the two species and the producer can invade effectively unimpeded by its competitor. We note that the minimum toxin strength required to produce an inhibition zone increases rapidly with increasing bacterial diffusivity, suggesting that even moderately motile bacteria must produce very strong toxins if they are to benefit in this way.

Figures

Tables

Table 1. Table 1: Eigenvalues σ 𝜎 \sigma and eigenvectors ( p ^ , s ^ ) ^ 𝑝 ^ 𝑠 (\hat{p},\hat{s}) of the Jacobian matrix ( 13 ), evaluated at the trivial and single-species stationary solutions ( p , s ) = ( p ∗ , s ∗ ) 𝑝 𝑠 superscript 𝑝 superscript 𝑠 (p,s)=(p^{*},s^{*}) of the classic Lotka-Volterra model given by the spatially-independent version of ( 7 )-( 8 ) with κ = 0 𝜅 0 \kappa=0 .

Table 2. Table 2: Eigenvalues σ 𝜎 \sigma and eigenvectors ( p ^ , s ^ , a ^ ) ^ 𝑝 ^ 𝑠 ^ 𝑎 (\hat{p},\hat{s},\hat{a}) of the Jacobian matrix ( 3 ), evaluated at the trivial and single-species stationary solutions ( p , s , a ) = ( p ∗ , s ∗ , a ∗ ) 𝑝 𝑠 𝑎 superscript 𝑝 superscript 𝑠 superscript 𝑎 (p,s,a)=(p^{*},s^{*},a^{*}) of the spatially homogeneous version of ( 7 )-( 9 ).

Equations

$$\frac{\partial P}{\partial T}$$

$$\frac{\partial S}{\partial T}$$

$$\frac{\partial A}{\partial T}$$

$$p=\frac{c_{PP}}{b_P}P,s=\frac{c_{SS}}{b_S}S,a=\frac{c_{PP}\delta A}{b_P\beta },$$

$$x=\sqrt{\frac{b_P}{D_A}}X,t=b_{P}T,$$

$$b$$

$$\kappa$$

$$ϵ$$

$$\frac{\partial p}{\partial t}$$

$$\frac{\partial s}{\partial t}$$

$$\frac{\partial a}{\partial t}$$

$$p^{\ast }=\frac{1-b{c}_s}{1-c_p{c}_s},s^{\ast }=\frac{b-c_p}{b(1-c_p{c}_s)}.$$

$$c_p<b<\frac{1}{c_s},$$

$$c_p>b>\frac{1}{c_s}.$$

$$J_{0}(p,s)=\mathopen{}\mathclose{{}\left(\small\begin{array}[]{cc}1-2p-bc_{s}s&-bc_{s}p\\ -c_{p}s&b(1-2s)-c_{p}p\par\end{array}}\right).$$

$$\sigma=-\frac{1}{2}\mathopen{}\mathclose{{}\left(p^{*}+bs^{*}\pm\sqrt{(p^{*}-bs^{*})^{2}+4bc_{p}c_{s}p^{*}s^{*}}}\right).$$

$$\begin{array}{rl}{p}_{\pm }& :=\frac{(b{c}_s+\zeta \kappa )a_{\pm }}{b{c}_s+\zeta \kappa a_{\pm }},\\ s_{\pm }& :=\frac{1-a_{\pm }}{b{c}_s+\zeta \kappa a_{\pm }},\\ a_{\pm }& :=\frac{B\pm \sqrt{B^2-4b(1-b{c}_s)\zeta {\kappa }^2}}{2\zeta {\kappa }^2},\end{array}$$

$$B:=b(1-c_{p}c_{s})+\mathopen{}\mathclose{{}\left((b-c_{p})\zeta-bc_{s}}\right)\kappa.$$

$$Q=\frac{(\sqrt{b(1-b{c}_s)}-\kappa \sqrt{\zeta }{)}^2}{b{c}_s+\zeta \kappa },$$

$$c_p$$

$$c_s$$

$$b-\kappa$$

$$\displaystyle\frac{1}{b}\mathopen{}\mathclose{{}\left(1-\frac{\zeta\kappa^{2}}{b}}\right)$$

$$c_p$$

$$\displaystyle\frac{1}{b}\mathopen{}\mathclose{{}\left(1-\frac{\zeta\kappa^{2}}{b}}\right)$$

$$c_p$$

$$c_s$$

$$b-\kappa$$

$$\displaystyle\frac{1}{b}\mathopen{}\mathclose{{}\left(1-\frac{\zeta\kappa^{2}}{b}}\right)$$

$$b-\kappa$$

$$\displaystyle\frac{1}{b}\mathopen{}\mathclose{{}\left(1-\frac{\zeta\kappa^{2}}{b}}\right)$$

$$J(p,s,a)=$$

$$\displaystyle\hskip 10.00002pt\mathopen{}\mathclose{{}\left(\scriptsize\begin{array}[]{ccc}1-2p-bc_{s}s&-bc_{s}p&0\\ -c_{p}s&b(1-2s)-c_{p}p-\kappa a&-\kappa s\\ \mu&-\mu\zeta\kappa a&-\mu(1+\zeta\kappa s)\end{array}}\right).$$

$$b-\kappa$$

$$0$$

$$\begin{array}{rl}\underset{x\to -\infty }{\lim }(p,s,a)& =(1,0,1),\\ \underset{x\to \infty }{\lim }(p,s,a)& =(0,1,0).\end{array}$$

$$p=\mathopen{}\mathclose{{}\left\{\begin{array}[]{ll}1,&\quad x\leq 0,\\ 0,&\quad x>0,\end{array}}\right.$$

Full text

1 Introduction

Bacteria perform many functions beneficial to life. For example, they are a key part of the digestive process in humans [1, 2, 3, 4] and other higher animals [5, 6, 7], form vital symbiotic relationships with plants [8, 9], including many food crops [10], and decompose waste organic matter [11, 12]. On the other hand, bacterial infections are also responsible for all manner of diseases in plants and animals, and therefore much effort is expended in controlling such pathogens. In the modern era, this has principally taken the form of antimicrobial substances such as antibiotics or disinfectants [13, 14]. However, with the advent of widespread antimicrobial resistance, this approach is becoming less and less viable [15, 16, 17, 18]. Moreover, broad and narrow spectrum antibiotics do not discriminate between helpful and harmful species, and so have the negative consequence of wiping out populations of beneficial species also [19]. It is therefore desirable to develop new methods of pathogen control which are capable of inhibiting harmful species while preserving the extant, beneficial, bacterial ecosystem [19]. For example, patients recently treated with antibiotics are vulnerable to potentially deadly infection by resistant Clostridium difficile as the destruction of the resident microflora clears the way for invasion. However, such infections have been successfully treated by re-establishing the microbiome via a faecal transplant from a healthy donor [20, 21, 22].

The microbiome plays a key role in preventing invasion of its host by other species, most pertinently pathogens [23]. A flourishing resident population that has filled the available ecological niches will, in theory, leave no room for an invader to establish within the community [24]. However, as well as competing indirectly via requisition of resources, many bacteria also compete directly by producing toxins which inhibit the growth of their rivals [25, 26]. If this inhibition is sufficiently effective at suppressing the target species, a toxin-producing invader may then be able to occupy the vacated niche in its stead [27, 28]. Conversely, a toxin-producing resident may be able to prevent an invader establishing merely by ensuring the environment is sufficiently unfavourable [23, 28, 29]. Thus a complete understanding of bacterial invasion and competition must take into account both the competition for resources familiar to ecological theory, and also the production of toxins. Understanding how these two competitive effects interact is the focus of the present work.

Mathematical modelling of such toxin-mediated bacterial competition has to date focused almost exclusively on growth in chemostats, as pioneered by Lenski and Hattingh [30]; a good introduction is provided in the review by Hsu and Waltmann [31]. Although the earliest work modelled toxin production via inhibitory pairwise reactions [30] (which is in fact simply classic Lotka-Volterra competition), later studies included a new dependent variable representing toxin concentration [32]. Many of the surveyed works concern the existence and stability of constant solutions to ordinary differential equation models, although limit cycles have been observed in the chemostat [33], and in a three-strain model comprising interactions between a toxin producer and both resistant and susceptible competitors, although that work did not include toxin concentration explicitly [34]. Further studies have included the analysis of rare mutations [35] and a rather detailed model introduced by Levin [36] which incorporates bacteria in different phases of growth and the recycling of organic material from dead bacteria being recycled, although in this case the toxin is externally applied rather than produced by bacteria.

A key feature missing from the current body of mathematical research in this area is the effect of spatial structure [37, 38]. It has been experimentally observed that spatial structure can alter the ability of a species to successfully invade an established population [28] and to establish a viable colony in the presence of competitors [39], when compared to the well-mixed case. This is because, when populations are well-mixed, competition is global and the weaker species cannot find a foothold. In contrast, when spatial structure is preserved, colonies grown from a random initial seeding of the domain may have time to establish themselves and are therefore potentially mature enough to resist invasion when rival colonies come into contact [39, 37]. Although one study [40] has looked at a reaction-diffusion model of a chemostat, they focused only on the stability of spatially homogeneous solutions and did not consider those exhibiting spatial structure. The present work will remedy this by considering competition in a spatially extended system; furthermore, rather than a chemostat we shall consider a closed system, it being representative of many bacterial ecosystems found in nature and the laboratory.

Microbial ecosystems are amazingly complex, incorporating not only bacteria but also protists, viruses, archea and fungi, and exhibiting all manner of antagonistic and mutualistic phenomena [41, 42]. When this ecosystem is located on or within a multicellular organism such as a human, then interactions between microbes and host add a further level of detail. Various microbes promote different host responses, promoting certain species and inhibiting others [43, 44, 23, 45, 46]. To understand these communities, we shall begin as simply as possible by considering a two-species ecosystem in a one-dimensional domain, with one species producing a toxin and the other susceptible to its effects.

After investigating the existence and stability of spatially homogeneous solutions corresponding to a well-mixed community, we allow the bacterial populations to vary spatially by adding one-dimensional diffusion. Inspired by the skin microbiome [43, 44], in which most species are nonmotile, we restrict our attention to weakly diffusing bacteria by considering the limit in which both bacterial diffusivities tend to zero. We analytically derive a steady solution to the spatial system in the case when the bacteria are completely immobile, i.e. when bacterial diffusivity vanishes, and derive conditions under which the toxin results in an inhibition zone near the producer within which the susceptible species cannot grow. Acknowledging that in reality even nonmotile bacteria exhibit some movement due to mechanical forces within the cell and in the local environment, we then solve our system numerically for small bacterial diffusivities. This enables us to describe invasion as a travelling wave connecting two single-species colonies, with one expanding its range as that of the other diminishes. We then interpret competitiveness in terms of the velocity of the travelling wave, and investigate how indirect competition for resources and direct inhibition through toxin production together determine a species ability to invade its competitor.

2 The producer-susceptible model

We consider two competing bacterial populations in a one-dimensional domain. One species, the ‘producer’, produces an antimicrobial toxin which inhibits the other, the ‘susceptible’. Denoting the concentrations of producer, susceptible and toxin at position $X$ and time $T$ by $P(X,T)$, $S(X,T)$ and $A(X,T)$, we model this ecosystem via the equations

[TABLE]

Here $b_{i}$ is the cellular division rate of species $i$ and $c_{ij}$ is the inhibitory effect of competition on species $j$ due to pairwise interactions with species $i$, for $i,j\in\{P,S\}$. The toxin is produced by $P$ at rate $\beta$ and degrades at rate $\delta$. $k$ is the inhibitory effect of the toxin on $S$, and $\eta$ the amount of toxin used up in inhibiting $S$. We assume that the producer, susceptible and toxin all move via random diffusion, with their respective diffusivity coefficients being $D_{P}$, $D_{S}$ and $D_{A}$. By using simple diffusion to represent bacterial movement, we are implicitly assuming that bacterial density is low enough that they move independently, as in higher density populations the ability of a bacterium to move may become restricted by its neighbours. Such a situation is better represented by biofilm models such as covered in [47] and is beyond the scope of this article. All coefficients in (1)-(3) are positive.

We nondimensionalise (1)-(3) by defining the dimensionless dependent variables

[TABLE]

dimensionless space and time

[TABLE]

and dimensionless constants

[TABLE]

Hence (1)-(3) is rendered

[TABLE]

This is the spatially extended Lotka-Volterra model [48, 49, 50] of two-species competition in one dimension, modified to explicitly account for the production of toxin by one species to inhibit the other. Note the somewhat unconventional nondimensionalisation (6); the more common choice in a Lotka-Volterra system is to write $(\hat{c}_{p},\hat{c}_{s})=(c_{p}/b,bc_{s})$. While this would allow us to take $b$ outside the brackets in (8) if we also wrote $\hat{\kappa}=\kappa/b$, the division rates then become conflated with the interspecific competition parameters. As the present work is focused primarily on competitive effects, we have chosen to keep these two processes separate for clarity.

3 Spatially homogeneous populations

We begin our analysis of (7)-(9) by seeking stationary and spatially homogenous solutions, i.e. by setting $p$, $s$, and $a$ to be constant. It is instructive to first consider the solution space with toxin strength $\kappa=0$, representing a completely immune target species $s$. In effect, this is equivalent to $p$ not producing a toxin at all, as $a$ is in this case immaterial to the population dynamics, and we can therefore neglect (9) and consider (7)-(8) only, without loss of generality. In the absence of spatial variation, this is the classic Lotka-Volterra model of two-species competition [51, 52]. There are up to four physically relevant (i.e. real and non-negative) stationary homogeneous solutions $(p,s)\equiv(p^{*},s^{*})$:

the trivial solution, $(p^{*},s^{*})=(0,0)$; 2. 2.

$p$ excludes $s$, $(p^{*},s^{*})=(1,0)$; 3. 3.

$s$ excludes $p$, $(p^{*},s^{*})=(0,1)$; 4. 4.

$p$ and $s$ coexist, with

[TABLE]

While the first three are always physically relevant, the coexistence solution requires either

[TABLE]

or

[TABLE]

The spatially homogeneous version of (7)-(8) with $\kappa=0$ has Jacobian matrix

[TABLE]

Eigenvalues $(\hat{p},\hat{s})$ and eigenvectors $\sigma$ for the trivial and single-species solutions are summarised in Table 1. We can therefore see that the trivial solution is unstable to perturbations in both $p$ and $s$; a population of $p$ only is always stable to perturbations in $p$ only, but is unstable to perturbations in $s$ provided $c_{p}<b$; a population of $s$ only is always stable to perturbations in $s$ only, but is unstable to perturbations in $p$ provided $c_{s}<1/b$.

The eigenvalues for the coexistence solution (10) are given by

[TABLE]

We omit the details of the eigenvectors for the sake of brevity; the salient detail is that each of the two eigenvectors has two non-zero components. Thus we see that $\sigma\in\mathbb{R}$ for all physically relevant parameter values, with one eigenvalue being positive only if $c_{p}c_{s}>1$. Hence the coexistence solution is stable if (11) holds, but not if (12) does. Combining this observation with the results of Table 1, we can now see that if (11) holds coexistence is the only stable solution, whereas if (12) holds then coexistence is an unstable manifold separating the two stable single-species solutions. If neither (11) nor (12) hold then the only stable stationary solution is single-species; producer only if $c_{p}>b$ and $c_{s}<1/b$, and susceptible only if $c_{p}<b$ and $c_{s}>1/b$. This is illustrated in Figure 1(a), depicting the solution space in the $(c_{s},c_{p})$-plane.

With the context of the classic system established, we extend the model to include toxin production by $p$ by setting $\kappa>0$ in the spatially independent version of (7)-(9). There now exist up to five physically relevant stationary solutions $(p,s,a)\equiv(p^{*},s^{*},a^{*})$:

the trivial solution, $(p^{*},s^{*},a^{*})=(0,0,0)$; 2. 2.

$p$ excludes $s$, $(p^{*},s^{*},a^{*})=(1,0,1)$; 3. 3.

$s$ excludes $p$, $(p^{*},s^{*},a^{*})=(0,1,0)$; 4. 4.

two possible branches of coexistence solutions, $(p^{*},s^{*},a^{*})=(p_{\pm},s_{\pm},a_{\pm})$, where

[TABLE]

with

[TABLE]

The first three solutions are always physically relevant, but the coexistence solutions (14) require $a_{\pm}\in\mathbb{R}$ and $0<a_{\pm}<1$. Defining

[TABLE]

and omitting the details for brevity, we find that $a_{+}$ is physically relevant provided

[TABLE]

or

[TABLE]

or

[TABLE]

and $a_{-}$ is physically relevant for either

[TABLE]

or

[TABLE]

Hence both solutions are physically relevant only for

[TABLE]

Note that the upper limit on $c_{p}$ in (15) increases monotonically from $b-\kappa$ to $b-\kappa/(1+\zeta\kappa)$ as $c_{s}$ increases across the given interval. Figure 1 illustrates these existence regions in the $(c_{s},c_{p})$-plane.

We can now analyse the stability of these solutions as for $\kappa=0$. The spatially homogeneous version of (7)-(9) has Jacobian matrix

[TABLE]

The eigenvalues $\sigma$ and eigenvectors $(\hat{p},\hat{s})$ of (3) evaluated at the trivial and single-species solutions of (7)-(9) are summarised in Table 2. Hence we see that the trivial solution is always unstable to perturbations in $p$ and $s$, but not to perturbations in $a$ only; a population of $p$ only is unstable to perturbations in $s$ provided $c_{p}<b-\kappa$, but is stable to perturbations in $p$ and $a$ only; a population in $s$ only is unstable to perturbations in $p$ provided $c_{s}<1/b$, but is stable to perturbations in $s$ and $a$. The trivial and single-species solutions of (7)-(9) therefore have similar existence and stability properties as when $\kappa=0$, the only difference being the shifting of the bifurcation line at which the $p$-only solution becomes stable to $c_{p}=b-\kappa$. Of course, if $\kappa>b$ then the $p$-only solution is stable for all positive values of $c_{p}$ and $c_{s}$, and there exists no region in the $(c_{p},c_{s})$-plane in which the $s$-only solution is the only stable solution; for this to occur in Lotka-Volterra dynamics, $s$ would require a negative birth rate.

There is a more striking qualitative difference between the solution space of (7)-(9) with $\kappa=0$ and that with $\kappa>0$, in that there is now a region (15) in which there are two possible coexistence solutions (14). While the stability properties of these solutions are analytically intractable, we can intuit them through consideration of the stability of the other solutions. For $c_{p}<b-\kappa$ and $c_{s}<1/b$, neither single-species population is stable; hence we expect the solution branch $(p_{-},s_{-},a_{-})$ to be stable here, as $(p_{+},s_{+},a_{+})$ is not physically relevant in this range. For $c_{p}>b-\kappa$ and $c_{s}>1/b$, it is $(p_{+},s_{+},a_{+})$ which is physically relevant, but not $(p_{-},s_{-},a_{-})$; here both single-species solutions are stable, so we expect $(p_{+},s_{+},a_{+})$ to be unstable, thus separating the two stable solutions. Finally, within the range (15) in which both coexistence branches are physically relevant, a population of $p$ is stable whereas a population of $s$ is not; we would therefore expect $(p_{+},s_{+},a_{+})$ to be unstable and $(p_{-},s_{-},a_{-})$ to be stable in this region also, with a saddle bifurcation along the line $c_{p}=b-\kappa+Q$, $(1-\zeta\kappa^{2}/b)/b<c_{s}<1/b$, on which the two solution branches intersect. This intuitive picture has been confirmed by numerical simulations for various parameter values in Matlab, the details of which we omit for brevity.

We can now see how the introduction of a toxin affects the classic picture of competition in the two-species Lotka-Volterra system. First, the bifurcation line at which a population of $p$ becomes stable, defined by $c_{p}=b-\kappa$, is moved closer to the $c_{s}$ axis as $\kappa$ increases from zero. In particular, for $\kappa>b$, the region of parameter space in which a population of $s$ is the only stable solution vanishes (see Figure (1)). Through production of a toxin, $p$ has enlarged the region of parameter space in which it can outcompete $s$. This effect is ameliorated somewhat by the level of inefficiency $\zeta$ of the toxin, where an increased value of $\zeta$ indicates a higher usage rate of the toxin. As $\zeta$ increases, the region of parameter space in which only a population of $p$ is stable decreases, being replaced by an expansion of the bistable region (15) in which both a population of $p$ and a coexisting population $(p_{-},s_{-},a_{-})$ are stable. In fact, as $\zeta\rightarrow\infty$, (15) becomes

[TABLE]

i.e. the region in which a population of $p$ is the only stable solution approaches $c_{p}>b$, $c_{s}<1/b$, coinciding with the equivalent region for $\kappa=0$. Note, however, the bifurcation line $c_{p}=b-\kappa$ is independent of $\zeta$. Therefore, production of a toxin always results in a solution space which is more favourable to $p$, as a larger region of the $(c_{s},c_{p})$-plane contains a stable population of $p$.

4 A two-colony solution with immotile bacteria

Now we have described the solution space of homogeneous solutions to (7)-(9), we turn our attention to solutions exhibiting spatial variation. As many bacteria are of low motility, for example those typically colonising human skin [43, 44] or otherwise adhering to host cells [53], we shall focus on the case when bacterial diffusivities are small and set $0\leq\epsilon\ll 1$. In the current section we shall allow the toxin to diffuse but set $\epsilon=0$. In this case we can find an exact stationary solution, which will aid our interpretation of travelling wave solutions to (7)-(9) with $0<\epsilon\ll 1$ in Section 5.

We choose to study solutions connecting a population of $p$ to a population of $s$ by imposing the boundary conditions

[TABLE]

Note that it does not matter at which end of the real line we choose to set the two different populations, as (7)-(9) are invariant under $x\rightarrow-x$. We shall see in Section 5 that, when $\epsilon>0$ and the bacteria are able to move, such boundary conditions allow us to cast an invasion process as a travelling wave connecting the two colonies. For now, with $\epsilon=0$, we can investigate how the production of a diffusing toxin by $p$ affects the growth of $s$ beyond the immediate vicinity of $p$, possibly resulting in an inhibition zone around the edge of the $p$-colony, within which $s$ is unable to grow.

With $\epsilon=0$, the steady versions of (7)-(8) become simple algebraic equations. We shall not consider coexistence solutions and instead focus on single-species populations; hence we have three possible solutions $(p,s)\in\{(0,0),(1,0),(0,1-\kappa a/b)\}$. We can therefore construct a solution by setting the bacterial concentrations to take one of these three options in successive regions of the real line, and then solving (9) to find the resultant toxin distribution. In order to satisfy the boundary conditions (20), we require a colony of $p$, $(p,s)=(1,0)$, in the leftmost region of the domain, and a colony of $s$, $(p,s)=(0,1-\kappa a/b)$, in the rightmost. The task at hand is therefore to find an appropriate solution connecting the two extrema which is consistent with equation (9). Although we can arbitrarily divide up the real line, for the sake of simplicity we shall assume that either the two colonies are in contact with one another, or are separated by a region in which both species vanish. We therefore have a colony of $p$, given by

[TABLE]

and a colony of $s$, given by

[TABLE]

for some $x_{0}\geq 0$. Note that we have exploited the translational symmetry of (7)-(9) to fix the edge of the $p$-colony, at which $p$ jumps from one to zero, to be at $x=0$. The edge of the $s$-colony, at which $s$ becomes non-zero, is located at $x=x_{0}$, where $x_{0}$ is an unknown to be found; if $x_{0}=0$, the colonies are in contact, whereas if $x_{0}>0$, they are separated by a region devoid of bacteria. We can therefore define the inhibition zone, in which $s$ cannot grow due to the inhibitory effect of the toxin, as the minimum possible value of $x_{0}$. Note that the inhibitory presence of the toxin diffusing from the $p$-colony results in $s$ being below its carrying capacity $s=1$ at the leftmost edge of the colony. Rather, $s$ tends to one as $x$ approaches infinity, as $a$ tends to zero in the same limit.

We now seek to determine $a$ by substituting (21) and (22) into (9), yielding

[TABLE]

This has solution

[TABLE]

where

[TABLE]

and we have exploited the requirement $a\rightarrow 0$ as $x\rightarrow\infty$ as per the boundary conditions (20) in order to derive the solution (26) in $x>x_{0}$. We choose the peak of the $\operatorname{sech}^{2}$-pulse in (26) to be at $x=x_{0}+x_{R}/m$ in order to simplify subsequent calculations.

It remains to determine the constants of integration $a_{L},a_{C},A_{C},x_{R}$ appearing in (26), and the location $x_{0}$ of the edge of the $s$-colony. This we achieve by requiring continuity of $a$ and its first derivative at $x=0$ and $x_{0}$, providing four equations in five unknowns. Thus one of the constants is arbitrary, albeit possibly subject to certain as yet undetermined constraints. This is a consequence of the immotile nature of the bacteria, as we shall see presently. Deferring the details to A for brevity, we choose $x_{R}$ to be our free parameter and hence arrive at the following expressions for the other four constants:

[TABLE]

Requiring that the solution is physically relevant, i.e. each of $p$, $s$ and $a$ are real and non-negative, then yields the upper bound on $x_{R}$

[TABLE]

where

[TABLE]

and $\hat{w}$ is the positive, real root of the cubic polynomial (A), arising from the requirement that $x_{0}\geq 0$. $\hat{w}$ is unique but exists only for toxin strength $\kappa<\kappa^{*}$, vanishing when $\kappa=\kappa^{*}$, and hence $x_{R}^{*}$ is piecewise continuous.

Now we have an upper bound on $x_{R}$, we can find the inhibition zone, i.e. the minimal value of the colony separation $x_{0}$, by minimising (28) subject to (29). As $x_{R}<0$, we find that (28) is always real, and has a minimum with respect to $x_{R}$ at

[TABLE]

But (31) is greater than the upper bound (29) on $x_{R}$, and so this minimum lies outside the physically relevant range of $x_{R}$. Hence $x_{0}$ is a monotonically decreasing function of $x_{R}$ for $x_{R}\leq x_{R}^{*}$. Therefore the inhibition zone is (28) evaluated at $x_{R}=x_{R}^{*}$, i.e.

[TABLE]

with the inhibition zone $I$ vanishing only if the toxin strength $\kappa\leq\kappa^{*}$, as defined in (30). Note that if $\kappa\leq 2b$ then $I=0$ for all $\zeta$.

Examples of two-colony solutions can be seen in Figure 2. In Figure 2(a) we see the colonies in contact, as the toxin strength is below the critical value $\kappa=\kappa^{*}$ and we have set $x_{R}=x_{R}^{*}$, its maximum. In Figure 2(b), we allow $x_{R}<x_{R}^{*}$, creating a gap between the two colonies. However, this is simply a consequence of our choice of the free parameter $x_{R}$ and is not in fact an inhibition zone. Because $\kappa$ is below the threshold, $s$ is capable of establishing itself closer to the producing colony, as in Figure 2(a), but is non-motile so cannot by itself move into the gap. A true inhibition zone is depicted in Figure 2(c), where the toxin strength has been increased past $\kappa^{*}$; here $x_{R}$ has again been set to its maximum, and therefore $x_{0}$ takes its minimum value, i.e. the size of the inhibition zone (32).

In Figure 3 we plot the inhibition zone $I$ against toxin strength $\kappa$ for various values of toxin inefficiency $\zeta$. Just above the threshold (30), we have

[TABLE]

so the inhibition zone at first grows linearly as $\kappa$ increases past $\kappa^{*}$, whereas for large $\kappa$

[TABLE]

increasing sublinearly. Hence the producer experiences diminishing returns for increasing toxin strength. Furthermore, an increase in either relative birth rate $b$ or toxin inefficiency $\zeta$ both increases the threshold (30) at which the inhibition zone becomes non-zero, and reduces the gradient $\,\mathrm{d}I/\,\mathrm{d}\kappa$ of the inhibition zone.

We conclude this section by considering the case $\zeta=0$, in which case the toxin concentration is depleted through natural degradation only, and not through interaction with $s$. In this case (23) still holds for $x\leq 0$, as the distribution of $p$ is unchanged, but (24) and (25) are identical equations. Therefore we have $a=A_{C}e^{-\sqrt{\mu}x}$ for all $x>0$; we no longer have the $\operatorname{sech}^{2}$ profile in (26) because it arises due to the nonlinearity in (25) provided by $\zeta>0$, and we set $a_{C}=0$ to prevent an unbounded solution as $x\rightarrow\infty$. Note that, since there is no longer a region in which $a$ is given by a $\operatorname{sech}^{2}$ profile, the solution family is in this case parameterised directly by $x_{0}$, rather than $x_{R}$. Imposing continuity of $a$ at $x=0$, and requiring non-negative solutions as above, we therefore obtain $a_{L}=A_{C}=1/2$, and

[TABLE]

and so the inhibition zone with $\zeta=0$ is given by the lower bound of (33), which is simply the limit of (32) as $\zeta\rightarrow 0$.

5 Travelling waves of weakly motile bacteria

The results of Section 4 hold for $\epsilon=0$, that is, bacteria which remain stationary for all time. This is of course an idealisation. In reality, processes such as cellular growth or interactions with the environment produce small forces which result in gradual movement even for those bacteria lacking specialised motility appendages such as flagella. Although it can be argued that the solution with $\epsilon=0$ is valid over short timescales, in order to gain a greater understanding of the bacterial dynamics we shall now consider $0<\epsilon\ll 1$, thereby allowing the bacteria to diffuse slowly throughout their environment. Motivated by the form of solution derived in Section 4 and illustrated in Figure 2, we shall continue to impose the boundary conditions (20). However, we now look for travelling wave solutions connecting a colony of $p$ at $x=-\infty$ to a colony of $s$ at $x=\infty$. Such a solution represents one possible invasion dynamic, in that we have an established population of each species, with both seeking to increase its colony size. A successful invader is thus one which outcompetes the other in the region where the two colonies meet, thus encroaching upon its rivals territory via diffusion; a useful background to such problems can be found in [54].

As is usual for a travelling wave analysis, we define the moving frame of reference $z=x-vt$, for some wavespeed $v$, and look for solutions to (7)-(9) which vary in $z$ only, yielding

[TABLE]

The boundary conditions (20) become

[TABLE]

Therefore, $p$ invading $s$ is represented by a wavefront moving to the right (with velocity $v>0$), and $s$ invading $p$ is represented by a wavefront moving to the left (with velocity $v<0$). By setting the boundary conditions (37) we have implicitly assumed that at least one of the single-species populations must be a stable solution of (7)-(9). Otherwise, the evolving dynamics will instead select the coexistence state $(p_{-},s_{-},a_{-})$ defined in (14), rather than a single species. Therefore we assume forthwith that at least one single-species solution is stable, i.e. that the competition parameters of (34)-(36) satisfy $c_{p}>b-\kappa$ or $c_{s}>1/b$ (see Table 2, the discussion around (3) and Figure 1). If there exists only one stable single-species solution, then clearly it is this which will invade the unstable population; we demonstrate this fact more rigorously in B by linearising the solution in the far-fields $z\rightarrow\pm\infty$ and showing that the velocity is bounded away from zero if precisely one single-species solution is stable, as per standard travelling wave theory [49, 55]. However, it is less clear what happens if both species exhibit stable populations, as in this case small perturbations in either far-field decay, and the velocity is therefore determined by fully nonlinear effects [56]. Indeed, this remains an active area of current research [57, 58], with recent progress made in the case of competitors with significantly differing diffusivities [59, 60, 61, 62]. This is in contrast to the current work in which we assume the bacteria have comparable motility. In order to make further progress, we shall now employ numerical methods to find travelling wave solutions to (7)-(9) with $\epsilon>0$.

We approximate the real line with far-field conditions (20) by a finite interval $x\in[0,d]$, for some $d>0$, with boundary conditions

[TABLE]

This approximation holds provided $d$ is sufficiently large that each dependent variable is exponentially close to constant-valued at the boundaries. We use second-order finite difference approximations to discretise the spatial derivatives in (7)-(9); the presence of the small parameter $\epsilon$ results in a stiff system, and so we integrate in time using the MATLAB function ode15s, optimised for such problems. Throughout we shall fix $b=D=1$, i.e. the two species have equal doubling times and diffusivities; we have carried out simulations with different values of $b$ and $D$, but we omit such cases as they do not significantly change the outcomes presented below. Instead, we focus our attention on varying the competition parameters $c_{p}$ and $c_{s}$, and the toxin parameters $\kappa$ and $\zeta$, in order to elucidate the relationship between classic Lotka-Volterra competition via pairwise interactions and direct inhibition through toxin production.

In each run of the simulation, we start from the initial conditions

[TABLE]

where $\Theta(x)$ is the Heaviside step function

[TABLE]

and $x_{1}\in[0,d]$. Most often $x_{1}=d/2$. The model is then integrated in time until a travelling wave solution is achieved; note this requires that $d$ is large enough that boundary effects are negligible throughout the simulation.

We plot examples of travelling wave behaviour for different values of $\epsilon$ and $\kappa$ in Figure 4, with the corresponding distribution profiles plotted in Figure 5. All four instances show the producer invading the susceptible even though $c_{s}>c_{p}$, with increases in $\epsilon$ or $\kappa$ corresponding to an increase in the speed of invasion; we see in Figure 4 the toxin diffusing ahead of the leading edge of its producer, clearing the way for expansion of the producing colony into space formerly occupied by the susceptible species. Of course, under different conditions the susceptible species may invade the producer. In Figure 5, we also plot the solution derived in Section 4 for $\epsilon=0$, taking the origin for this analytical solution to be the point where the simulated $p$ equals one half, in order to compare the profiles of bacteria with $\epsilon=0$ and $\epsilon>0$. In Figures 4(a)-(b) and 5(a)-(b), we see that there is no inhibition zone for $\epsilon=0$ (see (32)) and the solutions for $\epsilon=0$ and $\epsilon>0$ are qualitatively similar, with bacterial diffusion serving to smooth out the discontinuities exhibited when the bacteria are non-motile (see Section 4 and Figure 2). In contrast, in Figures 4(c)-(d) and 5(c)-(d), $\kappa$ is large enough that there does exist an inhibition zone for $\epsilon=0$. However, the size of this zone is greatly decreased even for $\epsilon=0.001$ (Figure 5(c)), and vanishes altogether for $\epsilon=0.01$ (Figure 5(d)).

In order to calculate the velocity of the moving frame of reference, we exploit the travelling wave framework (34)-(36). Integrating each of (34)-(36) across the real line, applying the boundary conditions (37) and changing the variable of integration from $z$ back to $x$ yields

[TABLE]

Starting from a particular choice of initial conditions, if the dynamics select a travelling wave solution we therefore have

[TABLE]

In the context of our numerical simulations, if $d$ is sufficiently large that $p$, $s$ and $a$ are asymptotically constant at the boundaries $x=0,d$, then (38)-(40) are good approximations of the velocity of the travelling wave. In practice, this requires integrating the system until each of the three expressions (38)-(40) are approximately equal and constant, checking that boundary effects remain negligible. We plot (38)-(40) for the parameter values of Figure 4(a) in Figure 6, showing each of the three expressions rapidly approaching an asymptote, the wavespeed.

To further elucidate the dependence of the velocity and inhibition zone on different parameters, we ran simulations of (7)-(9) over a range of values of the toxin strength $\kappa$, for nine different pairs of values of the competition parameters $c_{p}$ and $c_{s}$. The results can be seen for $\epsilon=0.001$ in Figure 7, and for $\epsilon=0.01$ in Figure 8. The presence of bacterial diffusion prevents us from defining the inhibition zone in the same way as in Section 4, as the bacterial distributions will now decay exponentially rather than vanish discontinuously. We instead define an effective inhibition zone as

[TABLE]

where $x_{\theta}$ is the region of space in which both $p$ and $s$ fall below the threshold $\theta$, i.e.

[TABLE]

We shall use $\theta=0.05$ throughout the present work.

In all cases shown in Figures 7 and 8, we see the velocity approaching an asymptote from below as $\kappa$ increases. Starting from $\kappa=0$, the inhibition zone is initially zero, but as $\kappa$ increases there comes a point at which it begins to grow. This point coincides with the magnitude of the velocity nearing its upper limit, with the inhibition zone thereafter increasing in size sublinearly as $\kappa$ increases. Thus we see that the velocity achieves its maximum value as the toxin becomes sufficiently inhibitory of $s$ near the leading colony edge of $p$ that an inhibition zone is formed. Once this has happened, $p$ is no longer actively competing with $s$ for access to resources, meaning that colony expansion depends only on the birth rate and diffusivity of $p$, rather than being impeded by the presence of $s$. In fact, if the inhibition zone is such that $s$ is exponentially small near the leading edge of the producing colony, $p$ is effectively governed by the Fisher-Kolmogorov equation [49]

[TABLE]

the leading-order approximation of (34) with $s\sim 0$. It is well known [49, 55] that the minimum wavespeed of (42) is $v=2\sqrt{\epsilon}$, and that the dynamics will always select the travelling wave profile associated with this value of $v$ for a wide range of physically relevant initial conditions. This formula gives $2\sqrt{0.001}=0.063$ (2 s.f.) and $2\sqrt{0.01}=0.2$, agreeing closely with the asymptotic values of $v$ plotted in Figures 7 and 8 as $\kappa\rightarrow\infty$. Thus we conclude that, once an inhibition zone has formed, the dynamics of $p$ are determined by the classic travelling-wave theory of (42). In contrast, although $s$ is also not directly competing with $p$ due to the inhibition zone, the toxin diffuses ahead of $p$ and impedes $s$, resulting in an effective growth rate of $b-\kappa a$ in (35) and allowing the producing colony to expand. Note the differing axes between Figures 7 and 8; this phenomenon requires much higher toxin strength as $\epsilon$ increases, due to the increased flux of susceptible bacteria towards the colony of producers.

For most data points in Figures 7 and 8, the velocity is positive, i.e. $p$ invades $s$. However, there are some exceptions in which toxin strength $\kappa$ is sufficiently low and $s$-competitiveness $c_{s}$ sufficiently high that $s$ is the invader. None of these instances correspond to $c_{s}=0.5$, as in this case the $s$ population is always unstable and so a travelling wave always moves so as to decrease its colony size. However, for $c_{s}>1/b$, where $b$ is the relative birth rate, the $s$ population is stable and so may be able to invade, provided it can outcompete $p$ and overcome the toxin. Note that this does not necessarily require instability of the $p$ population; in Figure 7(a) it can be seen that when $c_{s}=2.5$ there exists $\kappa>b-c_{p}=0.5$ such that $v<0$, i.e. both single-species solutions are stable but $s$ invades $p$. Of course, the majority of data points plotted over Figures 7 and 8 correspond to the bistable case with both single-species solutions being stable, but show $p$ invading $s$ due to the advantage of toxin production provided by our choice of $\kappa$ values. Note that there is no data point for $c_{p}=c_{s}=0.5$ and $\kappa=0$ in Figures 7(a) and 8(a), as in this instance both single-species solutions are unstable and the dynamics select the coexistence solution $(p_{-},s_{-},a_{-})$ defined in (14).

6 Discussion

Many bacterial species produce toxins as a means of increasing their fitness in highly competitive environments [25, 26]. By thus inhibiting its competitors, a toxin producer is afforded greater opportunity to exploit locally available resources. In order to understand this phenomenon in more detail, we have developed a model (1)-(3) of two species in competition. Note that, although the preceding analysis was carried out in terms of the nondimensional quantities (4)-(6), we shall in the present section discuss our results in terms of the original, dimensional quantities, in order to more effectively consider the biological reality.

Our model extends the classic Lotka-Volterra model in one dimension [48, 49] to explicitly include the production of a toxin by one of the two species. Indirect inhibition is modelled via the pairwise interaction terms $-c_{ij}ij$ in (1)-(2), where $i,j\in\{P,S\}$; if $i=j$ the term represents intraspecific competition, if $i\neq j$ it represents interspecific competition. Each competition constant $c_{ij}$ represents the inhibitory effect that the species $i$ has on the species $j$. Such effects depend not only upon species-specific traits such as rate of resource uptake, but also environmental factors such as precisely which resources are available, and how accessible they are. For example, if the two species depend primarily on different resources, then we would expect intraspecific competition to be greater than interspecific, and hence both $c_{ii}<c_{ij}$. On the other hand, a significant overlap of resource requirements with both species having approximately equal rates of uptake would result in each $c_{ij}$ to be approximately equal. If instead $P$, for example, had a greater rate of uptake than $S$, then we would have $c_{PS}>c_{PP}$ and $c_{SP}<c_{SS}$. Similarly, the birth rates $b_{i}$ will vary according to both species-specific and environmental factors. Therefore the location of any two-species system in parameter space depends not only on the choice of species, but also on the choice of environment. With this observation in mind, we shall now discuss our results in more detail.

Considering first spatially homogeneous solutions (i.e. a well mixed system), we showed in Section 3 that toxin production increases the region of parameter space in which a population of $P$ is stable. Moreover, if the inequality

[TABLE]

holds, then a population of $P$ is stable everywhere. Therefore toxin production increases the range of environmental conditions in which a species can dominate over its competitor, an especially useful trait when in a highly variable environment such as the skin or gut of a host animal. If the toxin strength $k=0$, then (43) can never hold, but even a very weak toxin opens up the possibility of the producer being stable in the whole of the physically relevant $(c_{s},c_{p})$-plane. Moreover, the inequality is independent of the interspecific competitiveness of $S$, meaning toxin production is effective against even highly competitive species, as the only means by which the susceptible species can negate (43) are increasing its birth rate $b_{S}$ or developing resistance to the toxin via decreasing its strength $k$. We note, however, that if (43) holds then there do exist bistable regions in which either the coexisting or $S$-only populations are also stable (see Figure 1(c)), and so the producer is not necessarily guaranteed to outcompete the susceptible even though it is stable solution.

The analysis of Section 3 provided a basis for the study of spatially varying solutions in Section 4, in which both the bacteria and the toxin diffuse across the domain. Motivated by skin bacteria such as Staphylococcus spp., we restricted our attention to weakly motile bacteria. In the special case of truly stationary bacteria, in which only the toxin diffuses, we derived an analytical, stationary solution describing adjacent colonies of the two species (see Figure 2). The producing colony acts as a toxin source, which diffuses into the susceptible colony, inhibiting its growth. Once the toxin strength increases above a certain threshold (30), this inhibition is total; an inhibition zone forms separating the two colonies. In dimensional form, (30) becomes

[TABLE]

i.e. the inhibition zone exists for $k>k^{*}$. Again, (44) is independent of interspecific competition. Note that (44) is quadratic in $b_{S}c_{PP}/b_{P}$; this nonlinearity arises from the depletion of the toxin as it interacts with the susceptible species, and is not present if the depletion coefficient $\eta=0$. Thus by decreasing the efficiency of the toxin, the susceptible species can increase the threshold at which the inhibition zone forms. This could be done either by increasing resistance, in that more toxin is required for the same inhibitory effect, or actively degrading the toxin molecules in the environment through secretion of some compound.

Upon allowing the bacteria to diffuse slowly in Section 5, we saw that two-colony solutions form travelling waves. Although the wave profiles are qualitatively similar to the discontinuous, stationary solutions of completely stationary bacteria—indeed, the stationary solution can be thought of as the singular limit of the travelling wave solutions as bacterial diffusivity $\epsilon\rightarrow 0$—two key differences were found. The first is that diffusion smooths out the discontinuities, so that bacteria distributions profiles decay exponentially rather than dropping to zero. The second, and more consequential, difference is that the threshold value of dimensionless toxin strength $\kappa$ required for an inhibition zone to form increases rapidly with bacterial diffusivity. Therefore, a more motile species must either invest much more in its toxin if it is to gain the benefits of an inhibition zone, or reserve expression of a toxin until it enters an immotile phase, such as pathogens adhering to host cells [53].

However, it is clear that there are substantial benefits to producing a toxin effective enough to form an inhibition zone. We saw in Figures 7 and 8 that the velocity asymptotes to a value independent of competition, and that the formation of the inhibition zone coincides with the velocity nearing its maximum magnitude. Thus the inhibition zone frees the producer from having to compete for the resources available at the leading edge of the colony; rather, the spreading speed is determined only by how effectively the producer can exploit those resources, and how motile it is. In such a scenario, the producer is then no longer under evolutionary pressure to compete with the inhibited species for resources, potentially allowing it to evolve increased competitiveness against other species in the community, or to become more efficient at exploiting secondary resources.

We note that we have not explored the possibility of travelling waves incorporating the coexistence solutions (14), in favour of focussing on the more tractable case of competing single-species colonies, the latter case being the more amenable to analysis. We defer study of the former to further work.

The present study provides insights that underpin our understanding of bacterial competition and presents a model of toxin-mediated inhibition to develop future investigations of species interactions, with potential implications for the therapeutic use of commensal bacteria [20, 21, 22, 63, 64].

Acknowledgements

The authors are grateful to the EPSRC for funding this work as part of the Liverpool Centre for Mathematics in Healthcare, under grant number EP/N014499/1. We are grateful to the anonymous reviewers for their constructive comments; in particular, for bringing to our attention the most recent literature on the theory of travelling waves in competitive systems.

Appendix A Calculation of the constants appearing in the two-colony solution of Section 4

In this appendix, we detail the derivation of the constants of integration $a_{L},a_{C},A_{C}$ and $x_{R}$ appearing in (26), and the location $x_{0}$ of the leftmost edge of the susceptible colony as defined in (22). This we do by imposing continuity of $a$ and its first derivative at the colony boundaries $x=0,x_{0}$. Requiring continuity of $a$ and its first derivative at $x=0$ yields

[TABLE]

with solution

[TABLE]

Then, requiring continuity of $a$ and its first derivative at $x=x_{0}$ yields

[TABLE]

We solve (46)-(47) for $a_{L}$ and $x_{0}$ in terms of $x_{R}$, obtaining

[TABLE]

and

[TABLE]

Thus (21)-(26) describe a family of stationary solutions to (7)-(9) with $\epsilon=0$, parameterised by $x_{R}$. We can now define the inhibition zone $I$ as the minimum value of $x_{0}$ with respect to $x_{R}$, i.e.

[TABLE]

It remains to find the range of values over which we can choose $x_{R}$, and then to determine the minimum of (49) within that range.

Requiring that $p$, $s$ and $a$ are all physically relevant imposes certain constraints that must be satisfied by our family of solutions. We note first that our solution (21) for $p$ is always physically relevant. Turning our attention next to $a$, requiring that (26) is physically relevant in $x\leq 0$ gives the condition

[TABLE]

For (26) to hold in $0<x\leq x_{0}$, with the constants of integration given by (45), we need

[TABLE]

as the left-hand side is minimal at $x=x_{0}$, we therefore require

[TABLE]

After rearrangement and comparison with (46), we see that (51) is always satisfied. Furthermore, if (51) holds then so does (50). Thus, because (26) is always real and non-negative in $x>x_{0}$, our solution for $a$ is physically relevant on the whole real line. However, from (22) we see that in order for $s$ to also be physically relevant we must have

[TABLE]

in $x>x_{0}$. The left-hand side of (52) is maximal at $x=x_{0}+x_{R}/m$, but the inequality fails at that point because $\alpha>b/\kappa$ (see (27)). Therefore the maximum of the left-hand side of (52) must be to the left of $x_{0}$, where $s\equiv 0$ and (52) does not apply. Hence we must have $x_{R}<0$, and the left-hand side of (52) decreases monotonically as $x$ increases from $x=x_{0}$. (52) is thus satisfied only if

[TABLE]

Finally, requiring $x_{0}$ as derived in (49) to be non-negative provides the condition

[TABLE]

Defining

[TABLE]

(53) indicates we must have $w\geq 0$; substituting for $\tanh x_{R}$ in (A) then gives

[TABLE]

By Descartes’ rule of signs, if the constant term in (A) is negative then the polynomial has precisely one positive root, which we denote by $\hat{w}$; otherwise it has no positive roots. As the inequality holds for large enough $w$, continuity requires that it also holds for $w\geq\max(0,\hat{w})$. Using (27) to determine conditions for negativity of the constant term in (A), for $x_{0}$ to be non-negative we must therefore have

[TABLE]

where

[TABLE]

and $\hat{w}$ is the positive, real root of (A), existing only for $\kappa<\kappa^{*}$ and vanishing when $\kappa=\kappa^{*}$.

Appendix B Some analytical results on the velocities of the travelling waves of Section 5

By recasting the system (7)-(9) in the moving frame of reference $z$, we have arrived at a system of ordinary differential equations, albeit with the unknown velocity $v$. We can make some analytical progress by linearising around the far-field solutions (37), following the method of Murray [49]. Considering first the right-hand limit $z\rightarrow\infty$, we write

[TABLE]

substitute into (34)-(36) and linearise in $e^{-\lambda_{R}z}$. Note that we fix $\lambda_{R}>0$ in order to ensure exponential decay as $z\rightarrow\infty$. After some minor rearrangement, this provides

[TABLE]

where

[TABLE]

Thus $v$ can be found in terms of $\lambda_{R}$ as an eigenvalue of (57), yielding three possible solutions for $v$ and the associated eigenvector $(\hat{p}_{R},\hat{s}_{R},\hat{a}_{R})$ for each choice of $\lambda_{R}$, namely

[TABLE]

where

[TABLE]

If the $s$ colony in the right-hand far-field is to connect smoothly to a non-zero value of $p$, it must be perturbed by a mode containing a non-zero $p$ component. The only eigenvector of (57) satisfying this requirement is (58), thus determining the velocity in terms of the perturbation growth rate $\lambda_{R}$. Although $\lambda_{R}$ is strictly positive by definition, its precise value depends on the choice of initial conditions. However, we note that if $c_{s}<1/b$, i.e. if a population of $s$ is unstable (see Section 3), then $v$ is positive for all positive values of $\lambda_{R}$. Therefore if a population of $s$ is unstable, a travelling wave solution satisfying (34)-(36) with far-field conditions (37) will always travel to the right, i.e. the producer will invade the susceptible. Furthermore, we can obtain a lower bound on the speed of the wave by minimising $v$ with respect to $\lambda_{R}$, namely

[TABLE]

We note this is equivalent to the linear spreading speed of van Sarloos [55].

A similar analysis can be performed for the $p$ colony in the left far-field. Letting $z\rightarrow-\infty$, we now write

[TABLE]

and linearise (34)-(36) as before. Thus we obtain

[TABLE]

where

[TABLE]

which we can solve to give the three possibilities

[TABLE]

where

[TABLE]

For the $p$ colony to connect to the $s$ colony, we now require a mode with non-zero $s$ component, the only candidate being (64). Hence we see that if $c_{p}<b-\kappa$ then $v$ is negative for all positive values of $\lambda_{L}$. In other words, if a population of $p$ is unstable then a travelling wave solution to (34)-(36) satisfying (37) will always travel to the left. Minimising $v$ with respect to $\lambda_{L}$ then gives an upper bound to the wavespeed of

[TABLE]

We therefore know that, whenever precisely one single-species solution is stable, a travelling wave connecting the two colonies must move so that the stable population invades the unstable one. Thus the strongly competitive species displaces the weakly competitive one. Note that (66) indicates that increasing $\kappa$ decreases the speed at which a weakly competitive producer is invaded; moreover, we saw in Section 3 that the region in which only the susceptible species is strongly competitive decreases in size as $\kappa$ increases, vanishing altogether when $\kappa=b$. With boundary conditions (37), we have a travelling wave moving to the right when only the producer is strongly competitive, and moving to the left when only the susceptible is strongly competitive; the question remains as to which direction the travelling moves when both species are strongly competitive. Comparing (58) and (64), we can again apply the logic of requiring all three components of both far-field perturbations to be non-zero to see that the velocity can only vanish when both single-species solutions are stable. Unfortunately, the above analysis is unable to elucidate any more than this simple fact, as the method requires the far-field solution to be unstable in order to derive a bound on the velocity.