Project 1.6. Plant-Herbivore System

Understanding

In the plant-herbivore problem, we examine the dynamics of introducing a population of 100 deer into an area of grassland and forest.
Our objective is to describe the types of dynamics that may arise (e.g. periodic populations, chaotic populations, etc.).

Plants will be measured in density units, whereas deer are to be modeled in terms of actual animals.

In the absence of deer, the plans observe a logistic growth with growth rate of .8 when population is low, and near 0 when the density is at 3000 density units. The initial plant population is at 3000 units.

The deer population declines at a rate of 1.1 when there is no food (i.e. loss is 1.1*pop(n) for the n+1 step). When food is present, there is population growth at a rate linearly increasing with plant density (0 at 0, and 1.5 when plants are at carrying capacity 3000 units). Of course, the plants are none too happy about this, and face a loss term linearly proportional to the deer population, 0 when there are 0 deer and increasing to 1.2 density units per deer when plants are at carrying capacity.

Plan

We will model the system using two discrete deterministic recurrence relations, which are coupled. We are asked to carry out three tasks:
Part I:
	1. Construct a simulation model and run it for 30 years under the following conditions
		a. Discrete time steps of one year
		b. Discrete time steps of one day
		c. Continuous time steps
	2. For each of the steps above, 
		a. construct time series of deer population and plant densities on the same graph.
		b. Also make a phase plane plot for each case (i.e. deer versus plants).
		c. Include comments about existence and stability of fixed points and differences between the various cases. 
	3. Which case is most  appropriate for this scenario?
	
Part II:
	Repeat Part I, only using the model of Caughley, in which the interaction terms are not linearly related. His model has decrease in plant density due to deer equal to 1.2D(1-exp(-.001*P)), and the increase in deer due to plants as 1.5D(1-exp(-.001*P)).

Part III:
	Compare the models in Part I and Part II. Plot (1-exp(-.001*P)) from 0 to 5000 and explain what this contributes to the second model over the first model. Which model does the best job of explaining "the observed behavior" (what "observed behavior"?).

Carry Out the Plan

Clear[P,De,n,R,K,Years];
R=.8;
K=3000.0;
Dnot=100.0;
Pnot=K;
Years=300;

P[0]=Pnot;
De[0]=Dnot;

P[n_]:=  P[n] =  P[n-1]*(R(1-P[n-1]/K)+1) - (1.2*P[n-1]/K)*De[n-1];
De[n_]:=De[n] = De[n-1] - 1.1*De[n-1] + (1.5*P[n-1]/K)*De[n-1];

Here's our simulation using the linear model of interaction:

plants=Table[{n,P[n]},{n,0,Years}];
deer=Table[{n,De[n]},{n,0,Years}];
both=Table[{P[n],De[n]},{n,0,Years}];
(* both//TableForm *)

Bonus step: thought that I'd present the growth models for the plant population for a variety of deer populations.

\(\(deerpop = 100; 

    Plot[{x, 
							
        x*\((R \((1 - x/K)\) + 1)\)\  - \ \((1.2*x/K)\)*.25*deerpop, 

        							
        x*\((R \((1 - x/K)\) + 1)\)\  - \ \((1.2*x/K)\)*.5*deerpop, 

        							
        x*\((R \((1 - x/K)\) + 1)\)\  - \ \((1.2*x/K)\)*deerpop, 

        							
        x*\((R \((1 - x/K)\) + 1)\)\  - \ \((1.2*x/K)\)*2*deerpop
		}, 

      	{x, 0, 3000}\ , 
	
      PlotStyle\  -> 
        \ {RGBColor[1, 0, 0], RGBColor[0, 1, 1], RGBColor[0, 1, 1], 
          RGBColor[0, 1, 1], RGBColor[0, 1, 1]}
	]; 

     (*\ 	PlotStyle\  -> 
        \ {Hue[{.2, .4, .6, .8, .9}], Thickness[.002], \ GrayLevel[0], 
          Dashing[{.01, .02, .03, .04, .05}]}
	*) \)\)

pplot=ListPlot[plants,PlotJoined->True,AxesLabel->{"Years",  "Density"},
	PlotRange->{1500,3000}];
dplot=ListPlot[deer,PlotJoined->True,AxesLabel->{"Years",  "Population"},
	PlotRange->{0,300}];
phaseplot=ListPlot[both,PlotJoined->True,AxesLabel->{"Plants",  "Deer"},
    PlotRange->{0,300}];

Show[pplot, dplot];

This time step (of a year) may be the most appropriate, since the deer are breeding once each year (not year round).
These are the plots of the various plant growth models under different conditions of deer population:

<<  Graphics`Graphics`\(DisplayTogether[
\ \ 
      Plot[x, \ {x, 0, 1}, PlotStyle\  -> \ {Hue[.6], PointSize[.01]}], 
\ \ 
      Plot[x*\((R \((1 - x/K)\) + 1)\)\  - \ \((1.2*x/K)\)*.25*deerpop, {x, 
          0, 3000}, 
\ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle\  -> \ {Hue[.2], PointSize[.03]}], 
\ \ 
      Plot[x*\((R \((1 - x/K)\) + 1)\)\  - \ \((1.2*x/K)\)*.5*deerpop, {x, 0, 
          3000}, 
\ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle\  -> \ {Hue[.4], PointSize[.05]}], 
\ \ 
      Plot[x*\((R \((1 - x/K)\) + 1)\)\  - \ \((1.2*x/K)\)*deerpop, {x, 0, 
          3000}, 
\ \ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle\  -> \ {Hue[.8], PointSize[.07]}], 
\ \ 
      Plot[x*\((R \((1 - x/K)\) + 1)\)\  - \ \((1.2*x/K)\)*2*deerpop, 
        \ {x, 0, 3000}, 
\ \ \ \ \ \ \ \ \ \ \ 
        PlotStyle\  -> \ {GrayLevel[0], Dashing[{.03}]}]
]\)

<< Graphics`MultipleListPlot`\(MultipleListPlot[plants, \ deer, PlotJoined\  -> \ {True, False}, 
	
      SymbolShape\  -> 

        \ \ \ {\ MakeSymbol[RegularPolygon[5, \ 3]], 
\ \ \ \ 
          MakeSymbol[

            \ \ \ \ \ \ {Line[{{2, \ 2}, \ {\(-2\), \ \(-2\)}}], 

              \ \ \ \ \ \ \ Line[{{\(-2\), \ 2}, \ {2, \ \(-2\)}}]}

            \ \ \ \ \ ]}
	]\)

Here's the simulation using the linear model of interaction, but daily intervals. Daily intervals has the effect of changing the growth rates, by dividing them up into 365 or so equal pieces.

Clear[P,De,plants,deer,both];
P[0]=Pnot;
De[0]=Dnot;
Delta=1/365.25;
Years=30/Delta;

P[n_]:=  P[n] =  P[n-1] + P[n-1]*R*(1-P[n-1]/K)*Delta - 
                    (1.2*P[n-1]/K)*De[n-1]*Delta;
De[n_]:=De[n] = De[n-1] - 1.1*De[n-1]*Delta + 
                    (1.5*P[n-1]/K)*De[n-1]*Delta;

plants=Table[{n,P[n]},{n,0,Years}];
deer=Table[{n,De[n]},{n,0,Years}];
both=Table[{P[n],De[n]},{n,0,Years}];
pplot=ListPlot[plants,PlotJoined->True,AxesLabel->{"Years",  "Density"},
	PlotRange->{1500,3000}];
dplot=ListPlot[deer,PlotJoined->True,AxesLabel->{"Years",  "Population"},
	PlotRange->{0,300}];
phaseplot=ListPlot[both,PlotJoined->True,AxesLabel->{"Plants",  "Deer"},
    PlotRange->{0,300}];

Here's the simulation using the numerical differentiation solver:

\(sol\  = \ 
      NDSolve[
\ \ \ \ \ \ {
          \(x'\)[t]\  == \ 
            x[t]*R \((1 - x[t]/K)\)\  - \ \((1.2  x[t]/K)\)*y[t], 
			
          \(y'\)[t]\  == \ \(-\ 1.1\) y[t]\  + \ \((1.5  x[t]/K)\)*y[t], 

          \ \ \ \ \ \ \ \ \ \ \ x[0]\  == 3000, \ \ y[0]\  == \ 100}, 
        \ {x, \ y}, \ {t, \ Years}]\)Plot[Evaluate[y[t]  /.  sol],  {t,  0, Years}]\(ParametricPlot[Evaluate[{x[t], \ y[t]}\  /. \ sol], 

      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {t, \ 0, \ Years}, \ 
      PlotRange\  -> \ All]\)

Here's the simulation using the exponential model of interaction:

Plot[1-Exp[-.001*x], {x,0,5000},PlotStyle -> {Hue[.6],PointSize[.01]}]

Clear[P,De,plants,deer,both];
P[0]=Pnot;
De[0]=Dnot;
P[n_]:=  P[n] =  P[n-1]*(R(1-P[n-1]/K)+1) - 
                   (1.2*(1-Exp[-.001*P[n-1]]))*De[n-1];
De[n_]:=De[n] = De[n-1] - 1.1*De[n-1] + 
                   (1.5*(1-Exp[-.001*P[n-1]]))*De[n-1];
Years=300;
plants=Table[{n,P[n]},{n,0,Years}];
deer=Table[{n,De[n]},{n,0,Years}];
both=Table[{P[n],De[n]},{n,0,Years}];
pplot=ListPlot[plants,PlotJoined->True,AxesLabel->{"Years",  "Density"},
	PlotRange->{400,3000}];
dplot=ListPlot[deer,PlotJoined->True,AxesLabel->{"Years",  "Population"},
	PlotRange->{0,1000}];
phaseplot=ListPlot[both,PlotJoined->True,AxesLabel->{"Plants",  "Deer"},
    PlotRange->{0,1000}];

Here's the exponential simulation using  daily intervals:

Clear[P,De,plants,deer,both];
Delta=1/365.25;
Years=30/Delta;

P[0]=Pnot;
De[0]=Dnot;

P[n_]:=  P[n] =  P[n-1] + P[n-1]*R(1-P[n-1]/K)*Delta - 
                   (1.2*(1-Exp[-.001*P[n-1]]))*De[n-1]*Delta;
De[n_]:=De[n] = De[n-1] - 1.1*De[n-1]*Delta + 
                   (1.5*(1-Exp[-.001*P[n-1]]))*De[n-1]*Delta;
                   
plants=Table[{n,P[n]},{n,0,Years}];
deer=Table[{n,De[n]},{n,0,Years}];
both=Table[{P[n],De[n]},{n,0,Years}];
pplot=ListPlot[plants,PlotJoined->True,AxesLabel->{"Years",  "Density"},
	PlotRange->{400,3000}];
dplot=ListPlot[deer,PlotJoined->True,AxesLabel->{"Years",  "Population"},
	PlotRange->{0,1000}];
phaseplot=ListPlot[both,PlotJoined->True,AxesLabel->{"Plants",  "Deer"},
    PlotRange->{0,1000}];

Here's the simulation using the numerical differentiation solver:

\(Years = 300; 

    sol\  = \ NDSolve[

        \ \ \ \ \ \ {
          \(x'\)[t]\  == \ 
            x[t]*R \((1 - x[t]/K)\)\  - \ 
              \((1.2*\((1 - Exp[\(-.001\)*x[t]])\))\)*y[t], 
			
          \(y'\)[t]\  == \ 
            \(-\ 1.1\) y[t]\  + \ 
              \((1.5 \((1 - Exp[\(-.001\)*x[t]])\))\)*y[t], 

          \ \ \ \ \ \ \ \ \ \ \ x[0]\  == 3000, \ \ y[0]\  == \ 100}, 
        \ {x, \ y}, \ {t, \ Years}]\)Plot[Evaluate[y[t]  /.  sol],  {t,  0, Years}]\(ParametricPlot[Evaluate[{x[t], \ y[t]}\  /. \ sol], 

      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {t, \ 0, \ Years}, \ 
      PlotRange\  -> \ All]\)

Stability of fixed point, linear case:

\(Clear[P, De, plants, deer, both]; 
P[0] = 1.1*K/1.5; 

    De[0] = R*\((K - P[0])\)/1.2; 
Years = 10; 


    P[n_] := \ \ 
      \(P[n]\  = \ \ 
        P[n - 1]*\((R \((1 - P[n - 1]/K)\) + 1)\)\  - \ 
          \((1.2*P[n - 1]/K)\)*De[n - 1]\); 

    De[n_] := 
      \(De[n]\  = \ 
        De[n - 1]\  - \ 1.1*De[n - 1]\  + \ \((1.5*P[n - 1]/K)\)*De[n - 1]\); 
    

Table[{P[n], De[n]}, {n, 0, Years}] // TableForm
\)\(Clear[P, De]; 
epsilon = .001; 
P[0] = 1.1*K/1.5*\((1 + epsilon)\); 

    De[0] = R*\((K - P[0])\)/1.2*\((1 + epsilon)\); 
Years = 300; 


    P[n_] := \ \ 
      \(P[n]\  = \ \ 
        P[n - 1]*\((R \((1 - P[n - 1]/K)\) + 1)\)\  - \ 
          \((1.2*P[n - 1]/K)\)*De[n - 1]\); 

    De[n_] := 
      \(De[n]\  = \ 
        De[n - 1]\  - \ 1.1*De[n - 1]\  + \ \((1.5*P[n - 1]/K)\)*De[n - 1]\); 
    

plants = Table[{n, P[n]}, {n, 0, Years}]; 

    deer = Table[{n, De[n]}, {n, 0, Years}]; 

    \(both = Table[{P[n], De[n]}, {n, 0, Years}];\)


     (*\ If\ you\ actually\ look\ at\ the\ values, \ 
      you\ see\ that\ the\ error\ is\ actually\ decreasing\ to\ 
        \(zero : \)\ *) 
 (*\ both // TableForm\ *) 

     (*\ let' s\ just\ look\ at\ the\ plots\ \(instead : \)\ *) 
\)\(pplot = 
      ListPlot[plants, PlotJoined -> True, 
        AxesLabel -> {"Years", \ \ "Density"}, 
	
        PlotRange -> {400, 3000}]; 

    dplot = ListPlot[deer, PlotJoined -> True, 
        AxesLabel -> {"Years", \ \ "Population"}, 
	
        PlotRange -> {0, 1000}]; 

    phaseplot = 
      ListPlot[both, PlotJoined -> True, 
        AxesLabel -> {"Plants", \ \ "Deer"}, 
\ \ \ \ 
        PlotRange -> {0, 1000}]; \)

Stability of fixed point, exponential case:

\(Clear[P, De, plants, deer, both]; 
Delta = 1/365.25; 

    Years = 1/Delta; 

P[0] = 1000.0*Log[15/4]; 

    De[0] = 1.5*R*P[0]*\((1 - P[0]/K)\)/\((1.1*1.2)\); 


    P[n_] := \ \ 
      \(P[n]\  = \ \ 
        P[n - 1]\  + \ P[n - 1]*R \((1 - P[n - 1]/K)\)*Delta\  - \ 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \((1.2*\((1 - Exp[\(-.001\)*P[n - 1]])\))\)*De[n - 1]*Delta\); 

    De[n_] := 
      \(De[n]\  = \ 
        De[n - 1]\  - \ 1.1*De[n - 1]*Delta\  + \ 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \((1.5*\((1 - Exp[\(-.001\)*P[n - 1]])\))\)*De[n - 1]*Delta\); 


    Table[{P[n], De[n]}, {n, 0, 10}] // TableForm
\)\(Clear[P, De]; 
Years = 35/Delta; 
epsilon = .001; 

    P[0] = 1000.0*Log[15/4]*\((1 + epsilon)\); 

    De[0] = 1.5*R*P[0]*\((1 - P[0]/K)\)/\((1.1*1.2)\)*\((1 + epsilon)\); 


    P[n_] := \ \ 
      \(P[n]\  = \ \ 
        P[n - 1]\  + \ P[n - 1]*R \((1 - P[n - 1]/K)\)*Delta\  - \ 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \((1.2*\((1 - Exp[\(-.001\)*P[n - 1]])\))\)*De[n - 1]*Delta\); 

    De[n_] := 
      \(De[n]\  = \ 
        De[n - 1]\  - \ 1.1*De[n - 1]*Delta\  + \ 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          \((1.5*\((1 - Exp[\(-.001\)*P[n - 1]])\))\)*De[n - 1]*Delta\); 


    plants = Table[{n, P[n]}, {n, 0, Years}]; 

    deer = Table[{n, De[n]}, {n, 0, Years}]; 

    both = Table[{P[n], De[n]}, {n, 0, Years}]; 
	

    Table[{P[n], De[n]}, {n, 0, 10}] // TableForm


     (*\ If\ you\ actually\ look\ at\ the\ values, \ 
      you\ see\ that\ the\ error\ is\ actually\ decreasing\ to\ 
        \(zero : \)\ *) 
 (*\ both // TableForm\ *) 

     (*\ let' s\ just\ look\ at\ the\ plots\ \(instead : \)\ *) 
\)\(pplot = 
      ListPlot[plants, PlotJoined -> True, 
        AxesLabel -> {"Years", \ \ "Density"}]; 

    dplot = ListPlot[deer, PlotJoined -> True, 
        AxesLabel -> {"Years", \ \ "Population"}]; 

    phaseplot = 
      ListPlot[both, PlotJoined -> True, 
        AxesLabel -> {"Plants", \ \ "Deer"}]; \)

Evaluate

One difference in the two models is the appearance of the stability of the equilibrium. In the case of the linear models, it seems that the system settles down to a stable fixed point.  In the case of the exponential model it is not clear that the solution tends toward a stable fixed point.

In the case of the exponential models on a yearly discrete scale, a limit cycle seems to get set up. On the other hand, for the exponential models with continuous scale, the system appears to be spirilling in on a stable fixed point.

As for which model does the best job of explaining the observed behavior, the exponential model does not meet the requirement that the interaction terms have the requisite slopes when the deer population is at its carrying capacities. The slope of the (1-exp(-.001P)) term is  exp(-.001P)*1000, which is approximately 0.0000497871 at 3000:

\(\(Exp[\(-.001\)*3000]/1000
		\)\)