Mathematical treatment of competition theory

by Prof. Matt Miller, Department of Mathematics, University of South Carolina
email: miller@math.sc.edu

Based upon RH MacArthur. 1972. Geographical Ecology. Harper & Row. pp 34-37.

> restart: gc():

Volterra proposed two logistically growing species had rates of increase that were determined by their own as well as their competitors' abundance.

> eq1a:=diff(X1(t),t)=(r1*X1/K1)*(K1 - X1 -alpha*X2);

> eq1b:=diff(X2(t),t) =(r2*X2/K2)*(K2 -X2 -beta*X1);

Alpha and beta are measures of the intensity of the interaction between the species (competition coefficients). When alpha and beta are zero, the species do not interact, and they grow as independent logistic populations. When alpha and beta are greater than zero, competitors reduce the equilibrium densities below the carrying capacities K1 and K2.

>

In order to determine whether the competitors will coexist or whether one will outcompete the other, we ask if each species can invade habitats fully occupied by the other. The world is a mosaic of patches of habitat, some empty, some occupied by one or more species, constantly changing as species invade new patches or go extinct in others. The probability is low that two species would simultaneously invade the same patches, so coexistence presumably arises as one species invades a patch previously occupied by another species. Presumably the occupant is near equilibrium (K), and the invader is rare when this process begins. In the context of our model, invasion is possible only if the invader's population growth rate dX1/dt > 0 when the occupant's density X2 = K2 (and dX2/dt > 0 when X1 = K1).

> coexist1:=rhs(eq1a)>0;

> K1/(r1*X1) * coexist1;

Let's now examine the condition where X1 is rare and the competitor is at carrying capacity ( the patch is occupied by the competitor):

> subs(X1=0,X2=K2,");

>

>

Warning, incomplete string;  use " to end the string

We now repeat the calculations for the other species: growth of species 2 must be positive inorder for it to invade:

> coexist2:=rhs(eq1b)>0;

> K2/(r2*X2) *coexist2;

> subs(X2=0,X1=K1,");

>

>

Warning, incomplete string;  use " to end the string

These two expressions, 0 < K1 - alpha K2, and 0 < K2 - beta K1, define the conditions under which each species can invade the other's habitat. They simplify to:

1/alpha > K1/K2 > beta

>

MacArthur suggests a vector equation solution to the same question:

At equilibrium

> eq3a:= K1 - X1 - alpha*X2 = 0;

> eq3b:= K2 - X2 - beta*X1 = 0;

> k1:=solve(eq3a,K1);

> k2:=solve(eq3b,K2);

To do matrix algebra we need to load the linear algebra package using this command:

> with (linalg):

Warning, the protected names norm and trace have been redefined and unprotected

Let's define two vectors for the interaction coefficients:

> vec1:=matrix(2,1,[1,beta]);

> vec2:=matrix(2,1,[alpha,1]);

Let's multiply each vector by the corresponding population density, and add them. Note that we have to evaluate the result of the scalar multiplication as a vector so we use the command :

evalm(vec1*X1) to show the result as a vector:

> evalm(vec1*X1) + evalm (vec2*X2);

Now we have two vectors that we are adding, so let's evaluate the sum as a vector:

> evalm(");

>

>

Warning, incomplete string;  use " to end the string

This looks a lot like the pair of expressions we wrote as k1 and k2 above:

> k1;

> k2;

Hence, the matrix expression is really the same as the pair of expressions for the equilibrium that we wrote earlier. Macarthur suggests that coexistence is possible only if the carrying capacities are a linear combination of the two interaction vectors (vec1 and vec2), in other words, the vector [K1,K2] must lie between the vectors [alpha,1] and [1,beta].

> with (plots):

Warning, the name changecoords has been redefined

Maple won't plot vectors directly so we have to convert them to lists after substituting in actual values for alpha and beta:

> vec1a:=convert(subs(alpha=0.8,beta=0.8,evalm(vec1)),vector);

> vec2a:=convert(subs(alpha=0.8,beta=0.8,evalm(vec2)),vector);

>

> v1:=plot([[0,0],convert(vec1a,list)],style=line,color=red):

> v2:=plot([[0,0],convert(vec2a,list)],style=line,color=green):

>

> display({v1,v2});

This graph is a phase plot, and the lines represent boundaries in phase space within which the carrying capacities must lie in order for the species to coexist.

>

We can plot the nullclines for this pair of species:

CHANGE THE PARAMETERS HERE IF YOU WANT TO TRY NEW VALUES FOR alpha AND beta OR K1 AND K2

> parms:=K1=1,K2=0.9,alpha=0.8,beta=0.8,r1=0.1,r2=0.15;

> nullcl_1:=solve(eq3a,X2);

> nullcl_2:=solve(eq3b,X2);

> nullcl_1a:=subs(parms,nullcl_1);

> nullcl_2a:=subs(parms,nullcl_2);

Here we plot the two nullclines. The $X1=0..1 command substitutes values from 0 to 1 in place of the variable X1 . Hence we get a plot of the points [0,(-K1-0)/alpha] and [1,(-K1 -1)/alpha] for nullcline 1.

> X1:=`X1`:pnc1:=plot([[X1,nullcl_1a] $X1=0..1],color=black):

> X1:=`X1`:pnc2:=plot([[X1,nullcl_2a] $X1=0..1],color=blue):

> kloc:=plot(subs(parms,[[K1,K2]]),style=point,color=black):

> vec1a:=convert(subs(parms,evalm(vec1)),vector):

> vec2a:=convert(subs(parms,evalm(vec2)),vector):

> v1:=plot([[0,0],convert(vec1a,list)],style=line,color=red):

> v2:=plot([[0,0],convert(vec2a,list)],style=line,color=green):

> display({v1,v2,pnc1,pnc2,kloc});

The red and green lines represent the boundaries of the carrying capacities within which coexistence is possible. The black dot shows the carrying capacities. The blue and black lines are the nulclines for species 1 and 2 respectively.

>

We'll now try a simulation with these parameters:

We will try several initial conditions:

X1=0.1,X2=0.1 (simultaneous invasion of empty habitat);

X1=K1,X2=0.1 (invasion of species 1 habitat by species 2);

X1=0.1,X2=K2 (invasion of species 2 habitat by species 1);

X1=K1,X2=K2 (instantaneous merging of two populations at K).

We use the subs(parms,[ ]) command to substitute our chosen parameter values into the list of equations and into the list of initial conditions.

> with(DEtools):

Warning, the name adjoint has been redefined

> dep:=DEplot(subs(parms,[eq1a,eq1b]),[X1(t),X2(t)],t=0..50,subs(parms,[[X1(0)=0.1,X2(0)=0.1],[X1(0)=K1,X2(0)=0.1],[X1(0)=0.1,X2(0)=K2],[X1(0)=K1,X2(0)=K2]]),scene=[X1,X2],linecolor=brown):

Let's plot the trajectories, vector field plot and nullclines on the same graph:

> display({dep,pnc1,pnc2});

>

What happens to the population trajectory when it crosses a nullcline?

How do the nullclines relate to the vector field plot?

Where is the equilibrium in relation to the nullclines and the population trajectories?

Is there more than one equilibrium?

What differences are there between equilibria associated with coexistence and equilibria associated with domination by one species? Look at the nullclines and the vector fields around the equilibria.

>

Try some different values for K1 and K2 or alpha and beta in the parms statemement above, and see how the nullclines and the range of possible K values change. What do the nullclines look like when there is coexistence? What do they look like when coexistence is impossible?

>