# Introduction to Mathematical Modeling in Ecology & Evolution

# Sally Otto

## Sponsors:
BIOS2 (NSERC - CREATE Training program)

Resources:
* An Introduction to Mathematical Modeling in Ecology and Evolution (2007) by Otto and Day
* Biomathematical modeling lecture notes (http://www.zoology.ubc.ca/~bio301/Bio301/Lectures.html)
* Mathematica labs (http://www.zoology.ubc.ca/biomath/labs.htm)
Maxima translation by Niki Love, Gil Henriques, and Sally Otto

# Part 1:  Classic one-variable models in ecology and evolution

## Exponential Growth 

### Assumes that each individual replicates at a constant rate over time:

 n[t + 1] = R n[t]           Discrete-time model ("Recursion equation")

dn/dt = n'[t] = r* n[t]     Continuous-time model ("Differential equation")

n[t]:   The population size at time t 
R:  The number of offspring per parent in a time unit (from t to t+1).
r:  The growth rate per parent per time unit.

## Logistic Growth 

### Assumes that each individual replicates at a rate that declines as a linear function of the current population size:

n[t+1] =(1+r (1-n[t]/K)) n[t]           Discrete-time model

dn/dt= n'[t][=r (1-n[t]/K) n[t]         Continuous-time model 

r:      The "intrinsic" growth rate per parent per unit time, realized when competition is weak (population size small)
K:      The "carrying capacity" defined as the population size at which competition causes the population to neither grow nor shrink

## Haploid Model of Selection

### Assumes that each individual carrying A has a fitness of 1+s relative to individuals carrying a, where p is the frequency of A:

 p[t+1] =((1+s)p[t]) / ((1+s)p[t] + (1-p[t]))       Discrete-time model

dp/dt = p'[t]=s*p[t] (1-p[t])                               Continuous-time model

p:  Frequency of allele A at a locus with two alleles (a and A).  p must lie between 0 and 1.

s:  Selection coefficient favoring allele A over a.

## Getting Used to Maxima

### Introduction to Simple Functions

Maxima can be used as a sophisticated calculator.

For example, if r = 0.3, K = 100, and n = 90 at some time, what population size do we expect in the next generation given logistic growth?
*Note: be sure to write kill(all) before each section, otherwise Maxima will remember parameter values from earlier sections. 

kill(all);

EQ1: (1 + r*(1-n/k))*n;

EQ1, r=.3, k=100, n=90;

### Plotting Functions in Maxima

For example, the per capita number of offspring per parent in the logistic model:

EQ2: (1 + r* (1 - n/k));

Plot this using the 'plot2d' function. Writing 'wx' before 'plot2d' keeps the graph in the Maxima window. You could plot by directly writing the equation with r=.3 & k=100:

wxplot2d([(1 + .3* (1 - n/100))], [n, 0, 250], [y, 0, 15], [xlabel, "n"], [ylabel, "kids"], [title, "Figure 1"]);

Or you could try a slightly more sophisticated method that will save time and help you keep track of parameters

EQ3:subst([r=.3, k=100],[EQ2]);

wxplot2d(EQ3, [n, 0, 250], [y, 0, 15], [xlabel, "n"], [ylabel, "kids"], [title, "Figure 2"]);

### Questions

The real power of Maxima, however, is that it can handle commands involving variables and parameters without having to specify their values.

For example, at what population size do we see the largest total number of offspring produced in the logistic model?

kill(all);

diff(((1 + r* (1 - n/K)) *n), n);

solve(%=0, n);

### Question 1: At what speed does an allele rise in frequency?

Using the continuous-time model, plot the rate of change  s p(1-p)as a function of p (between 0 and 1).  Choose whatever numerical value of s you wish.

### Question 2: When is evolutionary change fastest?

Again, using the continuous-time model, dp/dt=s p[t] (1-p[t]), what value of p maximizes the rate of allele frequency change?  Use D[ ] to find this result analytically, for any possible value of s. 

# Part 2: Equilibria and their stability

## Equilibria

*NOTE: we will use σ to denote the equilibrium population size as Maxima does not allow superscript in text. 

E.g., in the logistic model, if the population size at time t is n[t] = σ , it will remain at this value only if n[t+1] =σ.

This gives us an equation that we can solve for σ , by replacing all instances of n with σ in  n[t+1] = (1+r*(1-n[t]/K))*n[t].

solve( n = (1 + r*(1 - n/K)) *n, n);

Similarly, in continuous-time, we determine what value of σ causes dn/dt= r*(1-n[t]/K)*n[t] to equal zero:

solve(0=( r*(1-n/k)*n), n);

### Question 3: What are the equilibria for the haploid model of selection?

Find the equilibrium of p[t] (ρ) for both the discrete-time and continuous time models:

Discrete-time model
p[t+1]= (1+s) * p[t] / ((1+s)*p[t] + (1-p[t]))

Continuous-time model:
dp/dt = p'[t] = s*p[t]*(1-p[t])

First, use Maxima to solve for ρ analytically, then carry out the calculations by hand to confirm your answer.

## Stability

A system may or may not move toward a particular equilibrium.

A local stability analysis determines whether a system that starts near an equilibrium moves toward (stable) or away from it (unstable).

The idea:  We approximate the dynamics of a model near an equilibrium and then determine how it moves.  To do so, we use an incredibly useful approximation tool:  the Taylor Series.

C:\Users\Nicola Love\Pictures\taylor.series.png

taylor ( f(x[t]), x[t], σ, 10);

For example, Sin[x] can be rewritten as the following function around the point 1/2:

taylor(sin(x), x, .5, 4);

'ratprint: false ' suppresses warning messages about converting floating numbers to rational numbers. 

'ratsimp' cleans up the output

ratprint:false;

const: ratsimp(taylor(sin(x), x, .5, 0));

lin: ratsimp(taylor(sin(x), x, .5, 1));

quad: ratsimp(taylor(sin(x), x, .5, 2));

wxplot2d([cos(x),sin(x)],[x,-2*%pi,2*%pi],[y,-1,1],
[gnuplot_preamble,"set key bottom"]);

wxplot2d([sin(x), lin, const, quad], [x, 0, %pi], [y, 0, 1.5], [yx_ratio, .5], [legend, "sin(x)", "linear", "const.", "quad."], [style, lines, [points, 1], [points, 1], [points, 1]], [point_type, bullet]);

In a local stability analysis, we assume that we are near an equilibrium, and use the Taylor series to approximate the recursion or differential equation, f[n], to linear order:

ratsimp(taylor(f(n), n, σ, 1));

Describing the recursion equation as a function (f), where
 n[t+1] = f [n[t]], the Taylor series tells us that:

n [t+1] ~ f [σ] + ( n[t] - σ) * f '[σ] ( note that f '[σ] is the derivative of f with respect to n, evaluated at σ : df/dn | n=σ) 

n [t+1] ~ σ + ( n[t] - σ) * f '[σ]

 ( f [σ] = σ because the recursion started at an equilibrium that started at an equilibrium returns the equilibrium 

 ( n [t+1] - σ) ~ ( n [t] - σ) * f '[σ]  ( moved σ to the left)

so the distance to the equilibrium  ( n [t+1] - σ) changes by a factor λ = f '[σ] reach generation. 

Describing the differential equation as a function (f), where dn/dt= f [n [t] ], the Taylor series tells us that 

dn/dt ~ f [σ] + ( n [t] - σ) * f '[σ]
dn/dt ~ 0 + (n[t] - σ) * f '[σ]   
( f [σ] = 0 because there is no change at an equilibirum )


d ( n [t] - σ) / dt ~ (n[t] - σ) * f '[σ] 
( because d ( n [t] - σ) / dt = d(n[t])/dt - d (σ)/dt = d(n[t])/dt - 0

so the distance to the equilibrium changes at a rate r = f '[σ] 

For example, consider the logistic model in discrete time.  When will the system approach σ = 0, causing the system to go extinct?

derivative: diff((1 + r*(1 - n[t]/K))*n[t], n[t]);

λ: subst([n[t] = 0], [derivative]);

For λ to lie between -1 and +1, r must lie between -2 and 0.  [Technically, because 1+r measures the number of offspring per parent when competition is weak, r cannot fall below -1 or it becomes biologically meaningless.]

When will the system approach σ = K, causing the system to remain stably at carrying capacity?

λ: subst([n[t]=K], derivative);

This will be greater than one if:    r<0            σ=K is an unstable equilibrium
This will be between 0 and 1 if:   0<r<1        σ=K is a stable equilibrium
This will be between -1 and 0 if: 1<r<2         σ=K is a stable equilibrium with damped oscillations
This will be less than -1 if:            r>2             σ=K is an unstable equilibrium with damped oscillations

### Question 4: What determines stability of the equilibria for the haploid model of selection?

Under what conditions is ρ=0 stable? 

Under what conditions is ρ=1 stable? 


You can use either the discrete-time or the continuous-time model:

     p[t+1] =((1+s)p[t]) / ((1+s)p[t]+(1-p[t]))                  Discrete-time model
    
    dp/dt= p'[t]=s p[t] (1-p[t])                   Continuous-time model

Note that, in the discrete-time model, the fitness of A individuals is 1+s times greater than the fitness of a individuals, where 1+s must be positive (fitness can't be negative).

# Part 3: Beyond equilibria

## General Solutions

Some models can be solved generally, allowing us to predict the future state at any time in the future and how this depends on the parameters.

There are many methods for finding general solutions (see Chapters 6 & 9), including iterating recursions to deduce a general rule and using recipes to solve differential equations (e.g., separation of variables).

Maxima makes it easy to find general solutions for many of the simpler models. 

Recursion EQ: 

load("solve_rec")$
solve_rec(n[t+1] = R*n[t], n[t], n[0]=n0);

radcan(solve_rec(p[t+1] = (((1+s)*p[t])/((1+s)*p[t] + (1-p[t]))), p[t], p[0]=p0));

Differential EQ:

eqn: 'diff(n(t), t) = R*n(t);
atvalue(n(t), t=0, n0); 
sol: desolve(eqn, n(t)); 

### Question 5: Show that the logistic model in discrete time does not have a solution, but the model in continuous time does and then plot this solution.

Use: n[t+1] =(1+r (1-n[t]/K))n[t]       Discrete-time model

dn/dt = n '[t] = (r * (1- n[t]/K)) * n[t]   Continuous-time model

## Simulations

Some models, however, are too complicated to solve.  Some cannot be solved (e.g., the logistic in discrete time), while others may not be worth the time needed to obtain a general solution when all that is needed are solutions to specific cases.

kill(al);

n[t] := n[t-1] + n[t-1]*r*(1-n[t-1]/k);
n[0]: 10;
r: 1.5;
k: 100;
recursion: makelist(n[t], t, 0, 100);

wxplot2d([discrete,recursion], [style, [points, 1]], [xlabel, "time"], [ylabel, "population size"]);

For r values between 0 and 2, we see a stable equilibrium, as we would expect from the local stability analysis of Overscript[n, ^]=k.  For 2<r, we see oscillations, for 2.58<r, chaotic dynamics are observed, while for 3<r the system is driven extinct.

### Question 6: Simulating extinction

Add an If[ ] statement to the above simulation to return zero if ever the population size is predicted to be negative.

Use this simulation to explore what happens with r = 3.01 

## Numerical Solutions

Differential equations can also be solved numerically using Maxima's rk command,
which can work even for models that cannot be solved analytically.

n0: 10;
r: 1.5;
k: 100;
logistic: rk(r*(1 - n/k)*n,n,n0,[t,0,25,0.1])$

wxplot2d ([discrete, logistic])$

# Part 4:  Example of building a model from scratch

Scenario:  Let's model the number of species on a large landmass over long periods of evolutionary time, where the number of species is primarily influenced by speciation and extinction events within the landmass.  If extinction risk, d, is constant per species but speciation rate declines from an initial value, b, exponentially with the number of species already present (because of competition for overlapping resources/niches), then develop a model that might describe the number of species, n, on the landmass.

# Part 5:  Extending to models with more than one variable

Material to come.

# Part 6:  Another example of building a model from scratch

Scenario:  Consider two types of individuals that can help one another to raise a brood (e.g., one type might be better at detecting predators and the other type at collecting food).  Assume that neither of the two types would be able to grow on their own, with the number of offspring per parent, 1-r1 and 1-r2, respectively, being less than one.  Through cooperation, however, the per capita growth of each type rises linearly with the number of the other type, with slopes \[Rho]1 and \[Rho]2.

# Appendix: Quick reference to Maxima

## Getting help

<describe(command)>, or <? command> (with a space after the question mark): Prints to the screen the documentation of a particular command, e.g. <describe(integrate)> prints the documentation of the integrate() function. 
<describe(string, inexact)> or <?? inexact> (with a space after the question marks): Prints to the screen a list of items documented in the manual which contain "string" as part of their name.
<apropos("string")> returns a list of Maxima commands that have "string" appearing anywhere in them. For example, <apropos("exp")> retruns a list of functions including expand(), exp(), and ratexpand().

## Basic commands

Operations
<*> Times command (2*3; gives 6). 
<^> Power command (2^3; gives 8).
<!> Factorial (3!; gives 6).
Parentheses
<[ ]> Denotes a list, e.g. [2,3,4]. Can also denote the argument of a recurrence equation, e.g. n[3].
<( )> Has two functions. It can place variables together, e.g. (1+x)/(1-x) takes 1+x over 1-x. It is also used to denote arguments of functions, e.g. exp(3) is %e^3.
Grabbing output
<%o#> Grabs the previous output number #
<%> Grabs the previous output regardless of number. It is also used for some built-in symbols, e.g. %pi, %e, etc.
Other basic commands
<:> Defines a variable. For example: r : 10*x assigns the value 10*x to the symbol (variable) r. You can then perform operations on r,
e.g. 2*r returns 20*x. See also: section "Writing equations in Maxima".
<$> Omits the output of a particular operation. For example, a : 3*2; assigns the value 6 to the symbol a and returns the value 6,
whereas 3*2$ doesn't returns anything (but still performs the assignment).
<numer> Used to return numerical (instead of symbolic) values, e.g. 100*%pi, numer; returns 314.15926...
Substitutions
<f, object1 = object 2> Tells mathematica to substitute object1 with object2 in the function, e.g. 3*x^2, x=2*y+z; gives 3*(z+2y)^2
<subst(object1 = object2, f)> Same as the above

## Avoiding conflict with Maxima

Maxima tends to use lowercase letters for functions, so you can use initial upper case names for your functions and variables.
If you refer to previous outputs using %, it can be difficult to know exactly what your previous entry was. It is safer to assign a name to the output and then refer to this name later.
For example:

MyDerivative : diff(a*sin(b*x), x);
wxplot2d(subst([a = 1, b = 3], MyDerivative), [x, 0, %pi]);

## Functions and constants in Maxima (a small fraction!)

Absolute value, Logarithm, Square Root, Trigonometry
<abs(x)> Takes the absolute value of x
<log(x)> Takes the natural log of x. Maxima doesn't have a built-in function for base 10 logarithm or other bases. You can calculate the base-n logarithm of x as log(x)/log(n);
<sin(x)>, <cos(x)>, <tan(x)> Trigonometric functions
<arcsin(x)>, <arccos(x)>, <arctan(x)> Inverse trigonometric functions
<sqrt(x)> Square root
In-built symbols:
<%e> The exponential constant, 2.71838... %e^x can also be invoked using exp(x).
<%i> The square root of negative 1
<%pi> 3.14159...
<inf> Infinity

## Writing equations in Maxima

<x : y> Sets x to y immediately and from then on. Use kill(x) to unassign x, e.g. plot1 : wxplot2d(x^2, [x,0,10])
<'x : y> Does nothing until x is called, at which point x is assigned the value y
<x = y> Defines an equation
<f(x) := y> Defines a function f, e.g. f(x) := x^2
<f(x)> This gives the function evaluated at x, e.g. f(3), numer; gives 9 in the above example

## A list of helpful commands

Memory
<kill(symbol1)> Clears variable or function definitions
<kill(all)> Clears all variable or function definitions from memory
Algebra
<expand(expr)> - expands an expression, e.g. expand((1+x)^2) returns x^2 + 2*x + 1
<factor(polynomial)> - Self explanatory. E.g. factor(x^2+2*x+1) returns (x+1)^2
<collectterms(eqn, term)> - Collects parts of an equation involving a term and factors them separately. E.g. collectterms(a-b+a*x-2*b*x+a^2*x^2, x); returns a^2*x^2+(a-2*b)-b+a
<ratsimp(expr)> Simplifies a rational expression, e.g. a polynomial. E.g. ratsimp((x+2)*(x-2)); returns x^2-4. You can continue simplifying the expression until no further change occurs by calling fullratsimp() instead.
<radcan(expr)> - meaning "radical cancelation", it is useful in simplifying expressions containing logs, exponentials and radicals. E.g. radcan( (exp(x)-1)/(exp(x/2)+1)); simplifies an expression that ratsimp() cannot simplify.
<sum(f, i, imin, imax)> - Sums a function f which depends on i, from imin to imax. E.g. sum(i+1, i, 1, 4) returns 14
Calculus
<diff(f,x)> - Takes the partial derivative of f with respect to x. E.g. diff(x^2+y*log(x), x) returns 2*x + 1/2
<diff(f,x,n)> - Takes the nth derivative of f with respect to x. E.g. diff(x^2+log(x), x, 2) returns 2 - 1/(x^2)
<integrate(f,x)> - Finds the indefinite integral of f with respect to x. E.g. integrate(log(x),x) returns x*log(x)-x
<intgrate(f,x,xmin,xmax)> - Finds the definite integral from xmin to xmax. E.g. integrate(log(x),x,1,6) returns 6*log(6)-5
Solving a recursion equation for y as a function of x (in this example, y(t+1) = x*y(t)), SYMBOLICALLY, with initial value y(0) = y0:
load("solve_rec")$
solve_rec(y[t+1] = x*y[t], y[t], y[0]=y0);
Solving a differential equation for y as a function of x (in this example, y'(t) = x*y(t)), SYMBOLICALLY, with initial value y(0) = y0:
eqn: 'diff(y(t), t) = x*y(t);
atvalue(y(t), t = 0, y0);
desolve(eqn, y(t))
Solving a recursion equation NUMERICALLY:
***************************************************************************************
Solving a differential equation NUMERICALLY: rk uses the Runge-Kutta method
rk(t-x^2,x,1,[t,0,8,0.1])$
Solving an equation symbolically:
<solve(eqn, vars)>, e.g. solve(x^2 = 3, x)
Can also solve a system of equations, e.g. solve([x^2=3, y+3=x], [x, y])
Solving an equation numerically within a given interval. For example:
<find_root(x^2 = 7, x, 0, 100);> returns 2.6457...
Plotting:
<wxplot2d(f, [x,xmin,xmax])> Plots f versus x on the interval [xmin,xmax]
<wxplot2d([discrete, list1, list2], [style, points])> Plots points with x-coordinate given by list1, and y-coordinate given by list2
Others
<makelist(f, i, imin, imax)> Makes a table in list format of the function f with i values that run from imin to imax. E.g. makelist(i+3, i, 1, 4) returns [2,3,4,5]