Competition Models

by Prof. Matt Miller, Department of Mathematics, University of South Carolina

email: miller@math.sc.edu

This discussion is adapted from Robert H. MacArthur, 1972. Geographical Ecology. Harper & Row. pp. 33 - 58.

Volterra developed the basic theory of competition, using the logistic equation as his framework. His premise was that the `effective carrying capacity' of the environment is reduced by the presence of competitors for resources. Therefore we start from the logistic equation:

> restart: `~`:=``: with (plots): with(DEtools):

> logist:= diff(X1(t),t) = r1* X1 * (K1 - X1) / K1;

Warning, the name changecoords has been redefined

and introduce a term to describe the reduction in carrying capacity caused by the presence of members of a second species in the habitat:

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

X2 represents the population size of the competitors and alpha represents the per capita reduction in carrying capacity of species 1 caused by competition with species 2. Alpha is called the 'competition coefficient' of species 2 acting on species 1. Sometimes it is written with subscripts, i.e. alpha21, indicating the effect of species 2 on the growth of species 1.

A similar equation can be written for the second species, whose carrying capacity is reduced by the presence of the first:

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

Beta represents the per capita reduction in carrying capacity of species2 caused by competition with species 1. For similarity with equation 1, it is sometimes written as alpha with subscripts, i.e. alpha12, indicating the effect of species 1 on the growth of species 2.

>

These two equations when examined simultaneously, describe the dynamics of two species in competition with each other. There are no explicit functions that serve as solutions to these simultaneous equations. Hence a variety of methods have been used to examine their predictions. We can use variations on 'equilibrium analysis' to examine the behavior of the system.

A basic question we can ask is whether the two species can coexist or whether one drives the other to extinction. To examine this question, we first need to think about the biological conditions for coexistence. The real world is a mosaic of patches of habitat which are colonized and vacated regularly. New habitat opens up after environmental disturbances caused by winter cold, summer heat, hurricanes, fires, glaciers, volcanos, locust swarms, herds of goats, humans and other destructive forces. It is rare that two or more species will colonize an empty patch simultaneously. Instead, usually one species will enter a patch, establish itself, and then later be subject to invasion by other species. Therefore for COEXISTENCE to be possible, invaders must be able to enter a fully occupied patch of habitat.

>

Mathematically, this means that the invader (for instance species 2) must have a positive growth rate (dX2/dt > 0) during the first stages of invasion when it is rare (X2 is a small number). Since it is invading a habitat fully occupied by another species (for instance species 1), we need for dX2/dt > 0 even when the population X1 is at or near the carrying capacity K1. We write this as:

> invade2:= rhs( pop2 ) > 0;

To manipulate this inequality Maple has to be told some things that you know in your sleep: for instance, that the populations, the growth rates, and the carrying capacities are positive quantities. WARNING! Maple puts an annoying little flag (~) on all variables that come with assumptions; try to ignore it. (The funny little chicken scratches up on the first line can prevent it from showing up, but you won't always remember that little trick.)

> assume( X1 > 0 ); assume( K1 > 0 ); assume( r1 > 0 );

> assume( X2 > 0 ); assume( K2 > 0 ); assume( r2 > 0 );

When species 1 is invading, species 2 is at, or near, its carrying capacity and species 1 (the invader) is rare. We seek to have dX1/dt > 0 in this case.

> invade1:= rhs( pop1) > 0;

Since the world presumably consists of patches of species 1, patches of species 2, and empty space, coexistence (or patches of both) can only occur if BOTH SPECIES CAN INVADE THE OTHER. So we are interested in values of alpha and beta that make both invasion conditions true. Let us rearrange to put alpha and beta on one side of their respective inequalities.

> condition2:= solve( {invade2}, beta); condition1:= solve( {invade1}, alpha);

The following command gets rid of the ugly braces (technically we are telling Maple to select the first operand of the expression, which just so happens to be a set that contains exactly one operand, namely the inequality).

> condition2:= op( 1, condition2);

> condition1:= op( 1, condition1);

Recall that these conditions were supposed to hold under the most extreme conditions: one population (essentially) zero, and the other (essentially) at carrying capacity. Strictly speaking, we should compute limits here, but substituting is easier.

> condition2:= subs( {X2 = 0, X1 = K1}, condition2 );

> condition1:= subs( {X1 = 0, X2 = K2}, condition1 );

Note that the right hand sides of the two inequalities are reciprocals of one another. Let's call K1/K2 by a new name, R, and rewrite the condtions that allow for mutual invasion.

> condition1:= subs( K1 = R * K2 , condition1);

> condition2:= subs( K1 = R * K2 , condition2 );

Let's think about the consequences of this. What if alpha = beta = 0.9?

> restriction:= subs( {alpha=0.9, beta=0.9}, {condition1, condition2} );

We can tell Maple to solve the second of these inequalities and put it together with the first one, but this turns out to be harder than you might think. Anyhow, at some point it really makes more sense to do the work by hand. It's pretty easy to see that the two species can coexist only if the ratio R of their carrying capacities lies between 0.9 and 1.1. In the real world, carrying capacities fluctuate dramatically, sometimes by orders of magnitude, because of weather, predators, etc. Therefore it is very unlikely that two species could remain that close in ratio of K's for very long. What if alpha = beta = 0.1?

> restriction:= subs( {alpha=0.1, beta=0.1}, {condition1, condition2} );

How about this set of conditions for coexistence? Given fluctuations in K, is coexistence more or less likely than the last example? Try some other values for alpha and beta, not necessarily with alpha = beta. What general pattern do you see about the likelihood of coexistence, and the magnitudes of alpha and beta? How did we originally define alpha and beta? How does this relate to your observations?

Let's now look at the range of conditions under which coexistence can occur. From the above inequalities, if alpha = beta, coexistence can happen if the ratio of carrying capacities is between alpha and 1/alpha. We can plot this over the range 0.1 < alpha < 1.0:

>

> plot( {alpha, 1 / alpha}, alpha = 0.1 .. 1.0 );

>

An alternative analysis of the competition equations uses equilibrium methods directly. In this case we set the population growth rate to zero for both species.

> pop1; pop2;

At equilibrium both of these rates are zero, so we can simplify by multiplying each equation by K/(r*X):

> p1:=rhs(pop1)*K1/(r1*X1); p2:=rhs(pop2)*K2/(r2*X2);

We examine a case where the K1 = 900, K2 = 800, and both alpha and beta are 0.5. To see on the phase plane where these equilibria lie, we need to substitute these values in to the appropriate equilibrium expressions. Notice that R = K1 / K2 = 1.125, which is safely between 0.5 = alpha and 2.0 = 1 / alpha. For species 1, the equilibrium is:

> pp1:= subs({K1=900, alpha=0.5}, p1);

In order to plot this on the phase plane (X2 vs. X1) we need to rearrange this so it is X2 as a function of X1, so we solve the above expression for X2:

> nullcline1:= solve(pp1, X2);

We look at the equilibrium for species 2 in terms of species 1:

> pp2:=subs({K2=800, beta=0.5}, p2);

> nullcline2:= solve(pp2, X2);

Now we can plot both of the nullclines. When the populations lie above and to the right of these lines, they grow, and when they lie below and to the left, they decline.

> regions:= plot({nullcline1, nullcline2}, X1=0 .. 1600, X2= 0 .. 1600, color=black): # save this graph for later!

> regions; # and look at it now!

>

Now let's do a numerical solution of the equations, using the same values for K1, K2, alpha, beta. We will pick r1=r2=0.1 for this attempt:

> rate_equations:= [subs({K1=900, alpha=0.5, r1=0.1}, pop1), subs({K2=800, beta=0.5, r2=0.1}, pop2) ];

> vars:=[X1(t), X2(t)];

> domain:= 0 .. 120;

> init1:=[0,20,50]; init2:=[0,1400,100]; init3:=[0,100,1400];

> PhasePortrait:=DEplot(rate_equations, vars, domain, {init1,init2,init3}, stepsize=1, scene=[X1, X2]): #save it

> PhasePortrait; # see it

>

Have a look at the trajectories as t -> infinity. Where do they lie relative to the zero growth lines on the previous graph? The following command might be revealing.

> display( { PhasePortrait, regions} ); # overlay the phase portrait and the regions defined by the nullclines

Try some other values for K1, K2, alpha, and beta and see if you can demonstrate that the invasion criterion we derived earlier is correct. Note that you will have to start with your invader being RARE, NOT ABSENT. When we did the invasion analysis we set the invader to zero; obviously you will need to start with a few individuals or you will never get any growth.

>

>

Next we try to account for alpha and beta in a more mechanistic. less magical, way. What one species does to the other to affect its growth rate may involve several different aspects. There may be overt competition, physical exclusion from prime feeding or nesting sites, for example. Or the competition may be covert, say in consumption of a common food resource. Let us consider a pair of consumer populations X1 and X2 which eat resources R1 and R2, depriving each other in the process. We begin with MacArthur's equation (3):

> pop1:= diff(X1(t), t)/X1 = C[1] * (a[11]*w[1]*R1 + a[12]*w[2]*R2 - T[1]);
pop2:= diff(X2(t), t)/X2 = C[2] * (a[21]*w[1]*R1 + a[22]*w[2]*R2 - T[2]);

X1 and X2 are the population sizes of the consumers

R1 and R2 are the population sizes of the two resources.

w1 and w2 are the weights of individuals of each of the two resources.

aij is the probability that a consumer of species i encounters and eats an individual of resource j.

T1 and T2 are the maintenance requirements of consumer species 1 and 2.

C1 and C2 represent the conversion of grams of prey into population growth.

In this model, the per capita rate of population growth is determined by the amount of resource that exceeds maintenance requirements.

We now need to consider the population dynamics of the resources. We assume that each grows logistically when alone, but has its growth reduced by its encounters with consumers. Thus growth rate = logistic - loss to predator 1 - loss to predator 2; see MacArthur's equation (4):

> res1:= diff(R1(t), t)/R1 = (r[1]/K[1]) * (K[1] - R1) - a[11] * X1 - a[21] * X2;
res2:= diff(R2(t), t)/R2 = (r[2]/K[2]) * (K[2] - R2) - a[12] * X1 - a[22] * X2;

Determine equilibria equations for the two resources:

> equil_R1:= rhs(res1) = 0; equil_R2:= rhs(res2) = 0;

Solve for the equilibrium population sizes of the two resources (this is equation (5) of MacArthur on p. 38):

> Req1:= solve(equil_R1, R1); Req2:= solve(equil_R2, R2);

Determine the equilibria equations for the consumers. Since C1 and C2 cannot be 0, we can divide through by these and simplify our equations.

> equil_X1:= rhs(pop1) / C[1] = 0; equil_X2:= rhs(pop2) / C[2] = 0;

Now we substitute the equilibrium values for the resource populations into the equations for the consumers. In Maple we use the subs command to replace all occurrences of R2 and R1 in the equil_X1 and equil_X2 equations with the values we obtained by solving for R2 and R1 in the resource equations. This gives us the equations near the bottom of p. 38 in MacArthur.

> eX1:= subs({R2=Req2, R1=Req1}, equil_X1); eX2:= subs({R2=Req2, R1=Req1}, equil_X2);

We collect the constant terms, and the terms including X1 and X2, as done at the bottom of p. 38 and the top of p. 39. To understand the Maple commands convert and parfrac go to the help browser and use the keyword search. Or enter the commands ?convert, ?parfrac. Since we are doing nothing more than rearraging the equations, we might as well keep the same names.

> eX1:= convert(eX1, parfrac, {X1, X2}); eX2:= convert(eX2, parfrac, {X2, X1});

Note that what we got in both cases was: (-mess_1) * X1 + (-mess_2) * X2 + (mess_3) = 0,

which looks a lot like (mess_3) - (mess_1) * X1 - (mess_2) * X2 = 0,

which could be converted to (mess_3 / mess_1) - X1 - (mess_2 / mess_1) * X2 = 0,

which plainly has the form K1 - X1 - alpha * X2 = 0,

if (mess_3 / mess_1) = K and (mess_2 / mess_1) = alpha.

So let's try to extract the 'mess_1' by asking Maple to find the coefficient of X1 on the left hand side of the equation:

> b1:= coeff( lhs(eX1), X1);

Let's divide our original equation by this constant and collect terms:

> Eq6a:= convert( eX1 / b1, parfrac, {X1, X2});

Let's do it again for the other species. Get the coefficient of X2, divide through by this coefficient, and collect terms.

> b2:= coeff( lhs(eX2), X2);

> Eq6b:= convert( eX2 / b2, parfrac, {X1, X2});

Now we have a definition for alpha and beta, which are the coefficients in X2 in the first equation and the coefficients in X1 in the second equation:

> alpha:= coeff(lhs(Eq6a), X2); beta:= coeff(lhs(Eq6b), X1);

Let's examine these coefficients to see what they may signify. Remember that the parameters a11, a12 are the probabilities of consumption of resources 1 and 2 by consumer 1. Parameters a21 and a22 are the probabilities of consumption of resources 1 and 2 by consumer 2.

First consider the case where consumer 1 only eats resource 1 and consumer 2 only eats resource 2. This means that a12 = 0 and a21 = 0:

> A:= subs({a[12]=0, a[21]=0}, alpha); B:= subs({a[12]=0, a[21]=0}, beta);

What does this mean about the relationship between consumption of resources and competition?

Now let's consider the case where each consumer has 50% of its diet comprised of each resource:

> A:= subs({a[12]=0.5, a[21]=0.5, a[11]=0.5, a[22]=0.5}, alpha);
B:= subs({a[12]=0.5, a[21]=0.5, a[11]=0.5, a[22]=0.5}, beta);

Do you still believe what you thought about the relation between diet and competition? Consider consumer 1 eating 75% of resource 1 and 25% of resource 2, while consumer 2 eats 25% of resource 1 and 75% of resource 2:

> A:= subs({a[11]=0.75, a[12]=0.25, a[21]=0.25, a[22]=0.75}, alpha);
B:= subs({a[11]=0.75, a[12]=0.25, a[21]=0.25, a[22]=0.75}, beta);

Notice that under these conditions we do not have enough information to determine alpha and beta. If we assume the carrying capacities, growth rates and body sizes of the resources are similar then:

> A:= subs({a[11]=0.75, a[12]=0.25, a[21]=0.25, a[22]=0.75, K[1]=K[2], w[1]=w[2], r[1]=r[2]}, alpha);
B:= subs({a[11]=0.75, a[12]=0.25, a[21]=0.25, a[22]=0.75, K[1]=K[2], w[1]=w[2], r[1]=r[2]}, beta);

On the other hand if the carrying capacities and growth rates of the resources are different, then the competition coefficients are different, even with equivalent consumption. Here we assume K1 = 2*K2 and r1 = 0.5 * r2.

> A:= subs({a[11]=0.75, a[12]=0.25, a[21]=0.25, a[22]=0.75, K[1]=2*K[2], w[1]=w[2], r[1]=0.5*r[2]}, alpha);
B:= subs({a[11]=0.75, a[12]=0.25, a[21]=0.25, a[22]=0.75, K[1]=2*K[2], w[1]=w[2], r[1]=0.5*r[2]}, beta);

What generalization can you make from what you have seen about the relationship between diet overlap and the intensity of competition?

Look at alpha and beta again:

> alpha:= alpha; beta:= beta;

Some artful substitutions can make these complicated expressions look a lot simpler, and will suggest a natural generalization from 2 resources to an arbitrary number N.

> substitute_set1:= { a[11]^2 * w[1] * K[1] / r[1] = u[11]^2, a[12]^2 * w[2] * K[2] / r[2] = u[12]^2, a[21]^2 * w[1] * K[1] / r[1] = u[21]^2, a[22]^2 * w[2] * K[2] / r[2] = u[22]^2 };

> A:= subs( substitute_set1, alpha) ; B:= subs(substitute_set1, beta ) ;

If you take the square root on each side of each of the above substitution equations, you will notice that each term in the numerators of alpha and beta can also be writeen in terms of the u_ij's.

> substitute_set2:= {a[11] * a[21] * w[1] * K[1] / r[1] = u[11] * u [21] , a[22] * a[12] * w[2] * K[2] / r[2] = u[12] * u[22] };

> A:=simplify( subs( substitute_set2, A ) ); B:= simplify( subs( substitute_set2, B) );

Pretty clearly then, if there are 2 consumers and N resources (potentially) in common, we get competition coefficients ALPHA and BETA in terms of numbers U_ij that reflect the amount of resource sharing between the consumers and among themselves. Can you interpret what each U_ij term is measuring?

> ALPHA:= Sum(U[1][j]*U[2][j], j = 1 .. N) / Sum(U[1][j]^2, j= 1 .. N);
BETA := Sum(U[2][j]*U[1][j], j = 1 .. N) / Sum(U[2][j]^2, j= 1 .. N);

>