A 3-species competition model for bio-control

Applied Mathematics and Computation 218 (2012) 9690–9698

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A 3-species competition model for bio-control Ebraheem O. Alzahrani 1, Fordyce A. Davidson ⇑, Niall Dodds Division of Mathematics, University of Dundee, Dundee, DD1 4HN, Scotland, UK

a r t i c l e

i n f o

a b s t r a c t In this paper we model how the interaction of two species that compete for a common resource can be controlled by the introduction of a third species. This third species in the model can be considered as a bio-control agent. The aim here is to investigate whether the bio-control agent can slow, stall or even reverse the advance of one species into the territory of the other. We present conditions under which this can happen. The key point is that these conditions are on the relative strength of the bio-control agent. Thus, in application, these conditions could be met by choosing or designing the control agent appropriately and do not require the alteration of the properties of the existing species. Ó 2012 Elsevier Inc. All rights reserved.

Keywords: Reaction–diffusion Bio-control Competition Bistable

1. Introduction In this paper we model how the interaction of two species that compete for a common resource can be controlled by the introduction of a third species. This third species in the model can be considered as a bio-control agent. The aim here is to find conditions under which this bio-control agent can alter the advance of one species into the territory of the other. In order to do this, we first briefly discuss a class of competition models for which the interactions are assumed to be of classic ‘‘Lotka–Volterra’’ type, precise details will be given below (see also e.g. [1–3]). These models have the following structure for the interaction of two species u and v :

ut ¼ D1 uxx þ f ðu; v Þ;

v t ¼ D2 v xx þ gðu; v Þ

ð1Þ

þ

for ðx; tÞ 2 R  R . Typically, it is supposed that asymptotic conditions

ðu; v Þð1; tÞ ¼ S ;

ðu; v Þð1; tÞ ¼ Sþ ;

t>0

ð2Þ

are satisfied, for some suitable constant, non-negative vectors S . Here, x represents space and t, time. The restriction to one space dimension is made, in part, for clarity of exposition: the results here are not restricted to this case. Moreover, if the spatial domain is narrow or if the advance of either species is approximately planar, then this restriction is a reasonable model framework. The term ‘‘species’’ here is applied in its loosest sense and can be taken to represent different species of animals, microbes or plants. The species are assumed to move in a random manner (modelled via a diffusion term). The interaction of the species occurs through the functions f and g. In general, this can be synergistic, of predator–prey type or competitive. Here, we consider the last of these, where competition between the species is for a common resource, which is not explicitly modelled. Therefore, we consider kinetics for which fv < 0 and g u < 0. An interesting class of such models are those that are bistable, i.e. both S and Sþ are (asymptotically linearly) stable solutions of the associated kinetic problem. Certain conditions on f and g ensure that system (1) is monotone (see [4]). Thus, as in the scalar case, it can be shown that travelling wave solutions ⇑ Corresponding author. 1

E-mail address: [email protected] (F.A. Davidson). Current address: Department of Mathematics, Faculty of Science, King Abdulaziz University, P.O. Box 80200, Jeddah 21589, Saudi Arabia.

0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.03.012

E.O. Alzahrani et al. / Applied Mathematics and Computation 218 (2012) 9690–9698

9691

(i.e. solutions with constant speed and form) exist and can travel to the right, left or indeed be stationary. These travelling waves represent sharp transition fronts between regions where either u or v dominates to the exclusion of the other. Under the conditions studied here, these exclusion states are both stable and hence no a priori assumption can be made regarding the direction in which the front moves. However, in [5,6] details of how the interplay of relative motility (the ratio D : D2 =D1 ) and relative competitive strength determine the direction of travelling wave solutions are presented. (A more general and extensive discussion on travelling waves in such systems can be found in e.g. [4,7].) As a straightforward example, suppose D ¼ 1 and f and g have the typical forms f ðu; v Þ ¼ uð1  u  a1 v Þ; gðu; v Þ ¼ v ð1  v  b1 uÞ (see e.g. [2]). Then (1) has the trivial steady state ð0; 0Þ, a co-existence state ðu ; v  Þ and two semi-trivial states ðu; v Þ ¼ ð1; 0Þ and ð0; 1Þ. If 1 < a1 ; b1 , then it can be shown that the semi-trivial states are both stable. Setting S ¼ ð0; 1Þ and Sþ ¼ ð1; 0Þ, we have the following result: if 1 < a1 < b1 then a component-wise monontone left travelling wave exists for (1)–(2); if 1 < b1 < a1 then this wave travels to the right; if 1 < b1 ¼ a1 then the wave is stationary. If D – 1, the result is more complex and is determined by a combination of the system parameters that prescribe the relative potential energies associated with the end states S , see [5,6] for full details. Typical solution profiles are shown in Fig. 1a, where a sharp transition front is observed, travelling to the right for the parameter values chosen here. Behind the front, the wave rapidly relaxes to the (stable) steady state ðu; v Þ ¼ ð0; 1Þ and ahead of the front, ðu; v Þ ¼ ð1; 0Þ. For ease of visualisation, in Fig. 1b, the progression of the front is shown as a function of time for the species v only. Away from the interaction zone, u  1  v , so to the left of the front, v  1; u  0 and to the right v  0; u  1. For the simulation shown in Fig. 1 and all other simulations, numerical integration was performed using the MATLAB function pdepe, which invokes a finite difference approximation in space and utilizes the method of lines with an adaptive time step integrator for the resulting system of ordinary differential equations. The usual checks were performed to ensure that reducing the number of grid points and/or system tolerances did not alter the solution structure. We alighted on using 201 grid points, with absolute and relative tolerances set at 109 . No-flux boundary conditions were imposed at x ¼ 1; 1. Graphical output was also performed using MATLAB. If we now consider the introduction of a third species, w, system (1) becomes

ut ¼ D1 uxx þ f ðu; v ; wÞ;

v t ¼ D2 v xx þ gðu; v ; wÞ; wt ¼ D3 wxx þ hðu; v ; wÞ

ð3Þ

for ðx; tÞ 2 R  Rþ and it is supposed again that

ðu; v ; wÞð1; tÞ ¼ S ;

ðu; v ; wÞð1; tÞ ¼ Sþ ;

t>0

ð4Þ

for suitable constant vectors S . Three-component systems of this type have been studied previously. Of particular relevance to the work here, [8–10] discuss the stability of interacting fronts under the conditions that one species diffuses much slower than the other. Also, in recent related work, Lee et al. [11] discuss a more complex model for the bio-control in a 3-species predator–prey system, where, as here, the objective is to find conditions under which invasion can be altered. Patterns of patchy spread in three-species interaction models forms the focus in [13]. Bio-control is modelled in a 2-species model for a particular application in [14], see also references therein.

(a)

(b)

Fig. 1. Typical travelling wave solutions of the two-species competition model (1). (a) the species are given by u (blue), v (red) at 3 different time points. (b) The moving interface between v  1 (red) and v  0 (blue). Here f ðu; v Þ ¼ uð1  u  a1 v Þ; gðu; v Þ ¼ v ð1  v  b1 uÞ with a1 ¼ 3; b1 ¼ 2 and D ¼ 104 . As predicted, right travelling waves are formed (v is the stronger competitor). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

9692

E.O. Alzahrani et al. / Applied Mathematics and Computation 218 (2012) 9690–9698

For definiteness, here we will consider the functional forms f ðu; v ; wÞ ¼ d1 uð1  u  a1 v  a2 wÞ, gðu; v ; wÞ ¼ d2 v ð1  v  b1 u  b2 wÞ and hðu; v ; wÞ ¼ d3 wð1  w  c1 u  c2 v Þ. However, qualitatively similar results can be generated for other component-wise monotone kinetics, provided the associated kinetic system is tristable as defined shortly (and indeed for a larger class of models that follow the general form discussed in [5,6]). Note that w is assumed to compete for the same resource as u and v (this is its action of control), is self-sustaining and in the absence of u and v would propagate throughout the domain. The relative competitive strengths are measured by the parameters a1;2 ; b1;2 and c1;2 all of which are assumed to be positive constants. The relative motilities are given by ratios of the diffusion coefficients D1;2;3 . In the simulations shown below, we assume that both the motility and the intrinsic growth rates of the species are approximately equal and hence set D1;2;3 ¼ D and d1;2;3 ¼ d and then set d ¼ 1 as a new scaling of time. 2. The 3-component system 2.1. Tristability Consider the following ODE system associated with (3):

du ¼ uð1  u  a1 v  a2 wÞ; dt dv ¼ v ð1  v  b1 u  b2 wÞ; dt dw ¼ wð1  w  c1 u  c2 v Þ: dt

ð5Þ

This system is often referred to as the May–Leonard model as it is attributed to the authors who provided an extensive discussion of the dynamics of (5) in [12]. A general discussion of the existence and stability of steady states is provided but the main focus of that paper is on oscillations. Here, we restrict our attention to a special case as now described. After straightforward but lengthy calculations, it can be shown that system (5) has eight steady states: the trivial solution, six semi-trivial solutions and one co-existence solution given by

ðu ; v  ; w Þ ¼ ð1; 0; 0Þ; ð0; 1; 0Þ; ð0; 0; 1Þ;     1  a1 1  b1 1  a2 1  c1 ; ;0 ; ; 0; 1  a1 b1 1  a1 b1 1  a2 c1 1  a2 c1   1  b2 1  c2 0; and ; 1  b2 c2 1  b2 c2 

1  a1  a2 þ a1 b2 þ a2 c2  b2 c2 1  b1  b2 þ a2 b1  a2 c1 þ b2 c1 ; ; 1  a1 b1  a2 c1  b2 c2 þ a1 b2 c1 þ a2 b1 c2 1  a1 b1  a2 c1  b2 c2 þ a1 b2 c1 þ a2 b1 c2  1  c1  c2  a1 b1 þ a1 c1 þ b1 c2 : 1  a1 b1  a2 c1  b2 c2 þ a1 b2 c1 þ a2 b1 c2

ð6Þ

We are principally interested in the first three of these steady states. The states ð1; 0; 0Þ; ð0; 1; 0Þ represent the necessary ‘‘all of one and none of the other’’ dichotomy for the incumbent species far from the interaction zone. The state ð0; 0; 1Þ represents incumbent species removal by the bio-control. As is standard, the linear (asymptotic) stability of these steady states is determined by the eigenvalues of the Jacobian matrix: Jðu ; v  ; w Þ

0

1  2u  a1 v   a2 w B ¼ @ b1 v  c1 w

a1 u

a2 u

1  2v   b1 u  b2 w

b2 v 

c2 w

1  2w  c1 u  c2 v 

1 C A:

Again following standard calculations, it can be shown that provided all the interaction parameters are greater than one, then ð0; 0; 0Þ is an unstable node and ð1; 0; 0Þ; ð0; 1; 0Þ and ð0; 0; 1Þ are all stable nodes. Under these conditions on the parameters, it can also be shown that the other semi-trivial steady states are positive but unstable. Finally, additional conditions on the system parameters ensure that the co-existence steady state is positive (and hence biologically relevant), but is then also unstable. From now on, we assume that the interaction parameters are all greater than one and hence refer to system (5) as being tristable. Note, these steady states represent uniform steady states of the spatially extended system (3), which in this case we will also refer to as tristable.

9693

E.O. Alzahrani et al. / Applied Mathematics and Computation 218 (2012) 9690–9698

2.2. Initial conditions and asymptotic solutions We now discuss the solutions of the spatially extended system (3) with f ; g and h as given above. The kinetics suggest that travelling wave solutions can be anticipated that join the stable semi-trival states. Relevant to the application here are waves that satisfy w ! 0 as x ! 1. Moreover, we are interested in interaction fronts that separate zones of single species dominance and so below we consider wave solutions that connect ðu; v ; wÞ ¼ ð0; 1; 0Þ at x ¼ 1 to ðu; v ; wÞ ¼ ð1; 0; 0Þ at x ¼ 1, i.e. S ¼ ð0; 1; 0Þ and Sþ ¼ ð1; 0; 0Þ in (4). Indeed by applying similar arguments to those presented in [5], it can be shown that a component-wise monotone travelling wave solution of (3) exists. We now consider appropriate initial conditions associated with (3) and discuss how these evolve to form moving front solutions. First, we suppose that the u  v interaction is underway and that the bio-control agent is introduced at the site of interaction. An appropriate initial configuration for the system is required to represent this case. A combination of exponential functions can be used as a good approximation to the established u  v pair similar to that illustrated in Fig. 1a. It is assumed that w is introduced in a compact region and at a uniform density and hence a step function is used to model this. Fig. 2a depicts typical initial data. Notice, that this initial data is centred on x ¼ 0. In the simulations discussed below, the initial data used had the identical form but was either centred at x ¼ 0 or shifted to be centred at x  0:5. This shifting allowed for the evolution of the interaction over a longer time frame without interference with the numerically imposed boundary. Notice that in the absence of u and v, the initial p data ffiffiffiffi for w would evolve as two Fisher-type travelling fronts moving to the left and right, respectively with constant speed 2 D. Between the fronts, w  1. We now discuss how the solutions of the full system evolve. We return at the end to discuss the effects of varying this initial configuration, denoted by w0 ðxÞ, but note that the results are in general not sensitive to small quantitative changes in the initial configuration. In particular, we can replace the initial data for u and v shown in Fig. 2 with output from a numerical integration of the u  v only system without altering the results detailed below. However, for better comparison between cases and to clarify the discussions to come on the dependence on the initial configuration of w0 , the closed form approximations will be used throughout. 3. Results Numerical integration of the full system (3) using the initial data shown in Fig. 2 was conducted for a range of parameter values. Distinct classes of response were observed as is now discussed. Unless otherwise stated, in all that follows, parameter values were chosen such that in the absence of w, the u  v interaction would result in a right travelling wave as illustrated in Fig. 1 (i.e. a1 > b1 ). Corresponding results are obtained for the u  v left travelling wave case. 3.1. Annihilation of w In this section we will show that after some transient behaviour, w may be annihilated and u and before. This occurs under a number of conditions as now detailed.

(a)

(b)

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 −1

−0.5

0 x

0.5

1

0 −1

−0.5

0 x

0.5

v go on to interact as

1

Fig. 2. Initial conditions for system (3). The species are given by uðx; 0Þ (blue), v ðx; 0Þ (red) and wðx; 0Þ (green) where (a) bio-control: uðx; 0Þ ¼ 0:5 expð10xÞsechð10xÞ and v ðx; 0Þ ¼ 0:5 expð10xÞsechð10xÞ and wðx; 0Þ ¼ 1 for x 2 ½0:2; 0:2 and zero otherwise. (b) bio-buffer: uðx; 0Þ ¼ 0:5 expð10ðx  0:75ÞÞsechð10ðx  0:75ÞÞ and v ðx; 0Þ ¼ 0:5 expð10ðx þ 0:75ÞÞsechð10ðx þ 0:75ÞÞ and wðx; 0Þ ¼ 1 for x 2 ½0:2; 0:2 and zero otherwise. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

9694

E.O. Alzahrani et al. / Applied Mathematics and Computation 218 (2012) 9690–9698

Case A1: w contracts from both sides This behaviour is obtained if

a2 < c1 and b2 < c2 : In this simple case, w is out-competed by both u and Case A2: transient change of invasion speed If

c2  b2 > a2  c1 > 0;

ð7Þ

v and is annihilated. After a brief transient, u and v interact as before. ð8Þ

then the 3-component interaction zone moves to the right: v out-competes w and w out-competes u. However, the leading w  u front moves slower than the trailing v  w front. Hence, again after a transient period, w vanishes and u and v continue to move to the right. Note however that this transient could be of considerable length, depending on the relative speeds of the leading and trailing fronts. Moreover, the speed of each of these fronts is independent of the front speed of the u  v only interaction. Thus this transient can represent a period where the invasion of v into the territory of u is significantly altered. An illustration of this behaviour is given in Fig. 3, where the transient slows the invasion. Case A3: transient reversal of invasion If

0 < b2  c2 < c1  a2 ;

ð9Þ

then initially, the 3-component interaction zone moves to the left. Recall that in the absence of w, parameters are chosen such that v invades u, i.e, the travelling wave moves to the right. Hence, the introduction of w induces a reversal in the invasion: w out-competes v and u out-competes w, However, the leading v  w front moves more slowly than the trailing w  u front and hence w is eventually annihilated and this behaviour is only transient. The u  v interaction then proceeds to the right as before, see Fig. 4. 3.2. Permanent changes in dynamics In the above three cases, the annihilation of w means that the long-term behaviour is essentially determined by the u  v only system. We now discuss cases where w persists and therefore permanently alters the system dynamics. Case P1: permanent change of invasion speed This case happens if both leading and rear waves travel to the right and is obtained if the following conditions are satisfied:

0 < c2  b2 6 a2  c1 :

ð10Þ

If the equality is satisfied, then a fixed width interaction zone is formed, which travels to the right. An example of this behaviour is shown in Fig. 5a. As in the transient case detailed above, the speed of the interaction front can be different from that of the u  v only interaction and in Fig. 5a the invasion speed is much reduced (compare with Fig. 2). If the inequality is strict, then the interaction zone broadens in time with the leading front moving away from the trailing edge, see Fig. 5b. Special cases occurs when either c2 ¼ b2 or a2 ¼ c1 or both equalities are satisfied. In these cases, either the v  w or w  u or both interactions form standing waves. Thus the zone of interaction remains fixed in space at either one or both sides and the original competitors u and v are pinned and kept separate. In Fig. 5c, the case where both fronts are stationary

Fig. 3. Transient slowing of invasion. A transient 3-component interaction zone moves to the right at a reduced speed before w is annihilated and u  v interact as before. Numerical solutions of system (3) with a1 ¼ 3; a2 ¼ 2:01; b1 ¼ 2; b2 ¼ 2; c1 ¼ 2 and c2 ¼ 2:2; D ¼ 104 for t ¼ ½0; 500. (a) w  1 (red) w  0 (blue). (b) The moving interface between v  1 (red) and v  0 (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

9695

E.O. Alzahrani et al. / Applied Mathematics and Computation 218 (2012) 9690–9698

Fig. 4. Transient reversal of invasion. A transient 3-component interaction zone moves to the left before w is annihilated and solutions of system ð3Þ. (a) w and (b) v. Colour bar and parameter values as Fig. 3 except, a2 ¼ 2:0; c1 ¼ 2:5 and b2 ¼ 2:5.

v

invades u. Numerical

Fig. 5. Permanent change of invasion speed. Numerical solutions of system ð3Þ for t ¼ ½0; 500. In each case the plot shows w with w  1 (red) w  0 (blue). To the left of this zone, v  1; u  0 and to the right v  0; u  1. (a) A 3-component interaction zone of fixed width moves with reduced speed a2 ¼ 2:2; b2 ¼ 2; c1 ¼ 2 and c2 ¼ 2:2; (b) the interaction zone broadens a2 ¼ 2:5; b2 ¼ 2; c1 ¼ 2 and c2 ¼ 2:2; (c) segregation a2 ¼ b2 ¼ c1 ¼ c2 ¼ 2. All other parameter values are as in Fig. 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

is shown. Note that these conditions are independent of a1 and b1 . Hence, the introduction of the bio-control agent allows for spatial segregation of u and v irrespective of their relative competitive strength. In the 2-component system, such segregation is only possible in the special case a1 ¼ b1 . Case P2: permanent invasion reversal This case occurs if both leading and rear waves travel to the left and is obtained if:

b2  c2 P c1  a2 > 0:

ð11Þ

In this case, a region separating u and v moves to the left and hence the invasion is reversed. The speed of invasion is controlled by the relative competitive strengths of w and v (alt. u). If the equality is satisfied, then the interaction zone is of constant width, see Fig. 6. As in the case P1, there is a special case where c1 ¼ a2 and the trailing front is fixed. Thus u is pinned whilst w pushes back v.

3.2.1. Annihilation of v and u It may be that both the original species are unwanted and hence the role of the bio-control is to eradicate both. This can be obtained if w is the dominant competitor, namely:

0 < b2  c2

and 0 < a2  c1 :

See Fig. 7 for example behaviour.

ð12Þ

9696

E.O. Alzahrani et al. / Applied Mathematics and Computation 218 (2012) 9690–9698

Fig. 6. Invasion reversal. Numerical solutions of system (3) for t ¼ ½0; 500. In each case the plot shows w with w  1 (red) w  0 (blue). To the left of this zone, v  1; u  0 and to the right v  0; u  1. (a) A 3-component interaction zone of fixed width moves to the left a2 ¼ 2; b2 ¼ 2:2; c1 ¼ 2:2 and c2 ¼ 2; (b) the interaction zone broadens a2 ¼ 2; b2 ¼ 2:5; c1 ¼ 2:2 and c2 ¼ 2. All other parameter values are as in Fig. 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Species removal. The bio-control agent w dominates and annihilates both u and v. Numerical solutions of system (3) for t ¼ ½0; 500. The plot shows w with w  1 (red) w  0 (blue). To the left of this zone, v  1; u  0 and to the right v  0; u  1. a2 ¼ 2:2; b2 ¼ 2:2; c1 ¼ 2 and c2 ¼ 2. All other parameter values are as in Fig. 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

3.3. Dependence on inoculation 3.3.1. Inoculation pressure We return to discuss the effects of varying the ‘‘inoculation pressure’’ of the bio-control agent w. In the model, this is represented by varying the height and width (and therefore total ‘‘mass’’) of w0 ðxÞ. Extensive numerical experiments reveal that there is a complex relationship between the width and height of w and the relative competitive strengths of u; v and w. For the case where in the absence of w; u and v form a standing wave (a1 ¼ b1 ), the results are summarised in Fig. 8. For fixed relative strengths, provided the height and width of w0 ðxÞ are sufficiently large, then w persists and acts to alter the dynamics of the u  v system. Reducing the height and width seems to provide the intuitive relationship that the total inoculum mass requires to be sufficiently large: the boundary between persistence and annihilation is given approximately by mass = height  width = constant. However, there appear to be minimum heights and widths (Hmin and Wmin, respectively) for persistence: irrespective of the magnitude of the other dimension, if either the height or width falls below these minima, then w is eventually annihilated. These minima are functions of the other parameters (indicated by P in the figure), in particular they are increasing functions of values of the v  w and u  w interaction parameters a2 ; b2 ; c1;2 , irrespective of the relative level of these values. Note that in general, the minimum height is less than one. One might anticipate that this minimum height exists: if the height falls below the lowest (unstable) co-existence value as given in (6), then the system is likely to relax towards a configuration where w  0. However, this value does not seem to define the Hmin. For example, setting all competition parameters values a1;2 ¼ b1;2 ¼ c1;2 ¼ 2, from (6) the co-existence steady state value for w is

9697

Height

E.O. Alzahrani et al. / Applied Mathematics and Computation 218 (2012) 9690–9698

1

Persistance

Annihilation Hmin(P)

Wmin(P)

Width

Fig. 8. Dependence on Inocculation Pressure. The bio-control agent w either persists or is annihilated depending on the width and height of the initial inoculuum, w0 ðxÞ and shown schematically in the figure. The minima Hmin and Wmin and functions of other system parameters denoted here by P.

w  0:2. However, in this case, w is annihilated for significantly higher values (w0 ðxÞ 6 0:77), irrespective of the support of w0 ðxÞ. Dependence on the support of w0 ðxÞ could again be predicted as follows: if the support is too narrow, then the combined values of u and v in the support of w0 provides an additive competitive effect of sufficient strength to annihilate w. However, the initial data used here sums exactly to one for all values of x. Hence, this is unlikely to be the dominant process here. For example, again with all competition rates set at two, even if w0 ðxÞ ¼ 1 in its support, then w only persists provided this support is greater than 0.016. This value is well below the width of the interaction zone between u and v, which is approximately 0.5. What precisely defines these minima remains an open question. 3.3.2. Bio-buffers Finally, we consider an alternative initial configuration. As detailed above, Fig. 2a represents the situation where the extant species are interacting and the bio-control agent w is introduced at the zone of interaction. We could also consider the situation where the species u and v are yet to interact and the third species is placed inbetween them to act as a ‘‘bio-buffer’’, as illustrated in Fig. 2b. In this case, the v and u fronts advance towards each other, each with speed given by the Fisher speed pffiffiffiffi 2 D. With w0 ðxÞ  1 in its support and provided there are zones where all species are approximately zero separating w0 and u and v, then the resultant interactions are precisely as described above. However, as detailed immediately above, there is a complex dependence on the inoculum pressure. Moreover, the gap between the species seems to be an other complicating factor in deciding whether w persists and is therefore effective or is annihilated and the species interaction proceeds as before. For a fixed support of w0 , a relationship qualitatively similar to that schematically given in Fig. 8 exists between the height of w0 ðxÞ and the width of the gap between the fronts (assumed to be same on each side). In the absence of u and v, any non-zero compactly supported initial configuration of w would grow to a maximum height = 1 as the fronts spread as Fisher waves. However, for data as illustrated in Fig. 2b, but with w0 ðxÞ < 1 in its support, a balance between the time taken for w to increase the speed of the advancing v and u fronts occurs. If the fronts reach w before it has reached the minimum sustainable height as discussed above, then it is annihilated. Again this result is dependent on the choice of parameter values. 4. Conclusion In this paper, we considered the effect of introducing a third species to a classic 2-component bistable competition model. The role of this third species is as a bio-control agent, whose action is through competition for the common resource. The aim was to study conditions under which the population dynamics of the existing species could be altered by this agent. The model predicts that control of the interaction is dependent on there being sufficient inoculation pressure of the third agent. However, given this is in place, it was shown that under conditions on the relative competitive strengths of the species, the speed of invasion can be controlled and indeed, invasion can be halted or even reversed. Moreover, this control can be exerted on species that are already interacting or are yet to interact (in this latter case, the third agent can be thought of as a bio-buffer). A key feature of the model is that these changes in population dynamics can be induced without altering the relative competitive strengths or motility of the existing species. Rather, these results can be obtained by setting the competitive strength of the bio-control agent only. In practice, it would seem reasonable to assume that it would be very difficult to affect the dynamic response of existing species, but that it may be possible to ‘‘design in’’ or select the properties of any control agent, given knowledge of the incumbent species. The model results are therefore in line with these practical considerations. Many questions remain open including a fuller understanding of the required conditions on inoculum pressure and the effects of varying motility. These questions invite further numerical and analytical study. References [1] R. Cantrell, C. Cosner, Spatial Ecology via Reaction-Diffusion Equations, John Wiley and Sons Ltd, New York, 2003.

9698

E.O. Alzahrani et al. / Applied Mathematics and Computation 218 (2012) 9690–9698

[2] J. Murray, Mathematical Biology, II: Spatial Models and Biomedical Applications 2, Springer-Verlag, Berlin, 2003. [3] A. Okubo, S. Levin, Diffusion and Ecological Problems: Modern Perspectives, Springer-Verlag, New York, 2001. [4] A. Volpert, V. Volpert, V. Volpert, Traveling Wave Solutions of Parabolic Systems: Translations of Mathematical Monographs, 140, American Mathematical Society, Providence, R.I., 1994. [5] E. Alzahrani, F. Davidson, N. Dodds, Travelling waves in near-degenerate bistable competition models, Mathematical Modelling of Natural Phenomena 5 (5) (2010) 13–35. [6] E. Alzahrani, F. Davidson, N. Dodds, Reversing invasion in bistable systems, Journal of Mathematical Biology (2011), http://dx.doi.org/10.1007/s00285011-0490-9. [7] S. Petrovskii, V. Volpert, Reaction–diffusion waves in biology, Physics of Life Reviews 6 (2009) 267–310. [8] H. Ikeda, Multiple travelling wave solutions of three-component systems with competition and diffusion, Methods and Applications of Analysis 8 (3) (2001) 479–496. [9] H. Ikeda, Travelling wave solutions of three-component systems with competition and diffusion, Mathematics Journal of Toyama University 24 (2001) 37–66. [10] H. Ikeda, Dynamics of weakly interacting front and back waves in three-component systems, Mathematics Journal of Toyama University 30 (2007) 1– 34. [11] J. Lee, T. Hillen, M. Lewis, Continuous traveling waves for prey-taxis, Bulletin of Mathematical Biology 70 (2008) 654–676. [12] R. May, W. Leonard, Nonlinear aspects of competition between three species, SIAM Applied Mathematics 29 (1975) 243–253. [13] A. Morozov, S. Ruan, B.-L. Li, Patterns of patchy spread in multi-species reaction–diffusion models, Ecological Complexity 5 (2008) 313–328. [14] S. Petrovskii, K. McKay, Biological invasion and biological control: a case study of the gypsy moth spread, Aspects of Applied Biology 104 (2010) 37–48.