CHAPTER 7

1

1a

Define the formula for delta[p] and then plot it against p:

w11[p_]:=1-s p^2;
w12[p_]:=1-2s p(1-p);
w22[p_]:=1-s(1-p)^2;
wbar[p_]:=w11[p] p^2 + 2w12[p] p(1-p) + w22[p](1-p)^2;

delta[p_]:=p(1-p)(p(w11[p]-w12[p])+(1-p)(w12[p]-w22[p]))/wbar[p];

Plot[delta[p]/.s->0.1,{p,0,1}]

1b

Define the function in Equation 1 and then iterate it:

f[p_]:=(w11[p] p^2 + w12[p] p(1-p))/wbar[p]/.s->0.1

ListPlot[NestList[f,.2,500],PlotRange->All];

2

Define the phenotype frequencies at birth given the gene frequencies:

redfreq[{p_,q_}]:=(p^2+2p q+2p r)/.r->1-p-q;
yellowfreq[{p_,q_}]:=(q^2+2q r)/.r->1-p-q;
greenfreq[{p_,q_}]:=r^2/.r->1-p-q;

Define the phenotype fitnesses given the gene frequencies:

wred[{p_,q_}]:=1-.1redfreq[{p,q}];
wyellow[{p_,q_}]:=1-.1yellowfreq[{p,q}];
wgreen[{p_,q_}]:=1-.1greenfreq[{p,q}];
wbar[{p_,q_}]:=wred[{p,q}] redfreq[{p,q}]+
				wyellow[{p,q}] yellowfreq[{p,q}]+
				wgreen[{p,q}] greenfreq[{p,q}];

Define the recursion function foir the gene frequencies and iterate it for 500 generations:

func[{p_,q_}]:={wred[{p,q}](p^2 + p q + p r)/wbar[{p,q}],
	(wred[{p,q}] p q + wyellow[{p,q}](q^2 + q r))/
		wbar[{p,q}]}/.r->1-p-q

ans=NestList[func,{.4,.2},500];

Plot the phenotype frequencies:

ListPlot[Map[redfreq,ans],PlotRange->All]

ListPlot[Map[yellowfreq,ans],PlotRange->All]

ListPlot[Map[greenfreq,ans],PlotRange->All]

3

Define the Jacobian matrix from Equations 7.8b and 7.10 and find its eigenvalues:

f[i_]:=s(p[1]^2 + p[2]^2 + p[3]^2 - p[i])p[i];
jac=Table[D[f[i],p[j]] - D[f[i],p[3]],{i,2},{j,2}];
eigen=Eigenvalues[jac];

eigen/.{p[1]->1/3,p[2]->1/3,p[3]->1/3}

4

4a

Define the recurrence relations and evaluate the absolute values of the eigenvalues:

s=0.1; t=0.3;

f[1][x_,y_]:=(1-s y)x/(1-s x y-s(1-x)(1-y))

f[2][x_,y_]:=(1-t(1-x))y/(1-t x(1-y)-t(1-x)y)+
				m(1-2y)

jac=Table[D[f[i][p[1],p[2]],p[j]],{i,2},{j,2}];

mutation={0,0.003,0.01};
equil={{0.5,0.5},{0,0},{0,1}};

ans=Table[Eigenvalues[jac/.
	{m->mutation[[i]],p[1]->equil[[j,1]],
		p[2]->equil[[j,2]]}],{i,3},{j,3}]

Abs[ans]

4b

Define the recurrence relations in a form suitable for using in NestList and iterate with m=0; you should then examine and/or plot the output in a suitable way. Simulations with other values of m can be done in the same way.

g[{x_,y_}]:={f[1][x,y],f[2][x,y]}/.m->0

ans=ListPlot[NestList[g,{.4,.6},1000]];

5

5a

Solve the equalities in Equation 20 for the reward function in Equation 21:

w[1,1]=2/(n[1,1]+2n[1,2]);
w[1,2]=4/(n[1,1]+2n[1,2]);
w[2,1]=4/(n[2,1]+2n[2,2]);
w[2,2]=8/(n[2,1]+2n[2,2]);

Solve[{w[1,1]==w[2,1],
	w[1,2]==w[2,2],
	n[1,1]+n[2,1]==60,n[1,2]+n[2,2]==30},
	{n[1,1],n[1,2],n[2,1],n[2,2]}]

5b

Try to do the same for the reward function in Equation 22. The "solution" is impossible.

w[1,1]=2/(n[1,1]+n[1,2]);
w[1,2]=2/(n[1,1]+n[1,2]);
w[2,1]=4/(n[2,1]+2n[2,2]);
w[2,2]=8/(n[2,1]+2n[2,2]);

Solve[{w[1,1]==w[2,1],
	w[1,2]==w[2,2],
	n[1,1]+n[2,1]==60,n[1,2]+n[2,2]==30},
	{n[1,1],n[1,2],n[2,1],n[2,2]}]

Find the solution with all type-2 individuals in patch 2 and see if it satisfies Equation 20. It does.

ans=Solve[{w[1,1]==w[2,1],
	n[1,1]+n[2,1]==60,n[1,2]==0,n[2,2]==30},
	{n[1,1],n[1,2],n[2,1],n[2,2]}]

Transpose[Table[w[i,j],{i,2},{j,2}]]/.Flatten[ans]

6

The procedure func[n] takes the matrix n of n[[i,j]] values and returns the matrix of values in the next generation. sumn is the vector of the total number of individuals in the 3 patches (Equation 24); avpred is the vector of the average type of predator in the 3 patches (Equation 26); m is the matrix of interference coefficients (Equation 25); w is the matrix of rewards (Equation 23); next is the unstandardized matrix n in the next generation, and the last two lines standardize this matrix so that there are 100 individuals of each type.

func[n_]:=(sumn=Table[Sum[n[[i,j]],{j,10}],{i,3}];
	avpred=Table[Sum[j n[[i,j]],{j,10}]/sumn[[i]],{i,3}];
	m=Table[0.1 avpred[[i]]/j,{i,3},{j,10}];
	w=Table[i sumn[[i]]^-m[[i,j]],{i,3},{j,10}];
	next=Table[w[[i,j]] n[[i,j]],{i,3},{j,10}];
	sum=Table[Sum[next[[i,j]],{i,3}],{j,10}];
	Table[100 next[[i,j]]/sum[[j]],{i,3},{j,10}])

The simulation starts with equal numbers of individuals of each type in each patch and runs for 500 generations. The numbers in the final generation and their fitnesses are output.

start=Table[100/3,{i,3},{j,10}]//N;

res=NestList[func,start,500];

Transpose[res[[500]]]//TableForm

Transpose[w]//TableForm

7

I define the vectors of stay times, mating success and fitness:

time={10,20,40,60,100,150,200,250,300};
success={0.066,0.122,0.222,0.320,0.515,0.758,1.0,1.25,1.49};
fitness=success/(time+tau);

I define a function to plot fitness against time for a suitable range of values and for a given value of tau :

func[t_]:=(l=.99Min[fitness]/.tau->t;
			u=1.01Max[fitness]/.tau->t;
ListPlot[Transpose[{time,fitness}]/.tau->t,AxesOrigin->{0,l},
		PlotRange->{{0,300},{l,u}}])

This function can now be used to plot fitness against time for different values of tau, for example:

func[4]

8

Use pencil and paper.

9

I define the male and female frequency distributions and calculate their mean and variance:

maledata={1,1,6,1,5,4,15,18,26,25,17,33,18,36,23,
	21,17,8,8,6,1,1,1,0,1,0,1,0,0,0,1};
femaledata={0,0,0,0,0,0,1,2,2,1,0,4,7,5,11,14,14,
	17,11,12,9,5,10,14,8,3,2,1,0,0,1};

malemean=Range[31].maledata/Apply[Plus,maledata]//N

femalemean=Range[31].femaledata/Apply[Plus,femaledata]//N

malevar=((Range[31]^2).maledata - maletot malemean^2)/
	Apply[Plus,maledata]

femalevar=((Range[31]^2).femaledata - 
			femaletot femalemean^2)/Apply[Plus,femaledata]

I now define the normal density function, f[t], with the same mean and variance as the female distribution; F[t] is the corresponding cumulative normal distribution function, using the fact that F[z] = 0.5(1+Erf[z/20.5].

mu=femalemean;
sigma=Sqrt[femalevar];
f[t_]:=Exp[-(t-mu)^2/(2sigma^2)]/(sigma Sqrt[2 N[Pi]])
F[t_]:=0.5(1+Erf[(t-mu)/(sigma Sqrt[2])])

I define the function F[t] + f[t]/l and plot it against t to establish the approximate root of Equation 36, and then use FindRoot to obtain an accurate solution:

Remove[func]
func[t_]:=F[t]+f[t]/0.137

Plot[func[t],{t,0,50}]

tstar=t/.FindRoot[func[t]==1,{t,17.5}]

I define m[t] from Equation 35, plot it against t, and find the mean and variance of this distribution by numerical integration. (The part of the distribution m[t] with t < 0 is negligible and can be ignored.)

m[t_]:=If[t<tstar,f[t]+(D[f[u],u]/.u->t)/0.137,0]

Plot[m[t],{t,0,25}]

mean=NIntegrate[t m[t],{t,0,tstar}]

variance=NIntegrate[(t-mean)^2 m[t],{t,0,tstar}]