(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2"; MacintoshStandardFontEncoding; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8, 24, "Times"; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times"; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times"; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times"; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times"; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times"; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5, 12, "Courier"; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier"; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier"; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, 10, "Geneva"; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = leftheader, inactive, L2, 12, "Times"; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; fontset = leftfooter, inactive, L2, 12, "Times"; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; paletteColors = 128; automaticGrouping; currentKernel; ] :[font = section; inactive; preserveAspect; startGroup] CHAPTER 7 :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 1 :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] 1a :[font = text; inactive; preserveAspect] Define the formula for delta[p] and then plot it against p: ;[s] 3:0,0;23,1;31,0;60,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] 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; :[font = input; preserveAspect] delta[p_]:=p(1-p)(p(w11[p]-w12[p])+(1-p)(w12[p]-w22[p]))/wbar[p]; :[font = input; preserveAspect; endGroup] Plot[delta[p]/.s->0.1,{p,0,1}] :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] 1b :[font = text; inactive; preserveAspect] Define the function in Equation 1 and then iterate it: :[font = input; preserveAspect] f[p_]:=(w11[p] p^2 + w12[p] p(1-p))/wbar[p]/.s->0.1 :[font = input; preserveAspect; endGroup; endGroup] ListPlot[NestList[f,.2,500],PlotRange->All]; :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 2 :[font = text; inactive; preserveAspect] Define the phenotype frequencies at birth given the gene frequencies: :[font = input; preserveAspect] 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; :[font = text; inactive; preserveAspect] Define the phenotype fitnesses given the gene frequencies: :[font = input; preserveAspect] 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}]; :[font = text; inactive; preserveAspect] Define the recursion function foir the gene frequencies and iterate it for 500 generations: :[font = input; preserveAspect] 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 :[font = input; preserveAspect] ans=NestList[func,{.4,.2},500]; :[font = text; inactive; preserveAspect] Plot the phenotype frequencies: :[font = input; preserveAspect] ListPlot[Map[redfreq,ans],PlotRange->All] :[font = input; preserveAspect] ListPlot[Map[yellowfreq,ans],PlotRange->All] :[font = input; preserveAspect; endGroup] ListPlot[Map[greenfreq,ans],PlotRange->All] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 3 :[font = text; inactive; preserveAspect] Define the Jacobian matrix from Equations 7.8b and 7.10 and find its eigenvalues: :[font = input; preserveAspect] 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]; :[font = input; preserveAspect; endGroup] eigen/.{p[1]->1/3,p[2]->1/3,p[3]->1/3} :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 4 :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] 4a :[font = text; inactive; preserveAspect] Define the recurrence relations and evaluate the absolute values of the eigenvalues: :[font = input; preserveAspect] s=0.1; t=0.3; :[font = input; preserveAspect] f[1][x_,y_]:=(1-s y)x/(1-s x y-s(1-x)(1-y)) :[font = input; preserveAspect] f[2][x_,y_]:=(1-t(1-x))y/(1-t x(1-y)-t(1-x)y)+ m(1-2y) :[font = input; preserveAspect] jac=Table[D[f[i][p[1],p[2]],p[j]],{i,2},{j,2}]; :[font = input; preserveAspect] mutation={0,0.003,0.01}; equil={{0.5,0.5},{0,0},{0,1}}; :[font = input; preserveAspect] ans=Table[Eigenvalues[jac/. {m->mutation[[i]],p[1]->equil[[j,1]], p[2]->equil[[j,2]]}],{i,3},{j,3}] :[font = input; preserveAspect; endGroup] Abs[ans] :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] 4b :[font = text; inactive; preserveAspect] 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. ;[s] 3:0,0;64,1;72,0;225,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] g[{x_,y_}]:={f[1][x,y],f[2][x,y]}/.m->0 :[font = input; preserveAspect; endGroup; endGroup] ans=ListPlot[NestList[g,{.4,.6},1000]]; :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 5 :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] 5a :[font = text; inactive; preserveAspect] Solve the equalities in Equation 20 for the reward function in Equation 21: :[font = input; preserveAspect] 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]); :[font = input; preserveAspect; endGroup] 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]}] :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] 5b :[font = text; inactive; preserveAspect] Try to do the same for the reward function in Equation 22. The "solution" is impossible. :[font = input; preserveAspect] 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]); :[font = input; preserveAspect] 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]}] :[font = text; inactive; preserveAspect] Find the solution with all type-2 individuals in patch 2 and see if it satisfies Equation 20. It does. :[font = input; preserveAspect] 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]}] :[font = input; preserveAspect; endGroup; endGroup] Transpose[Table[w[i,j],{i,2},{j,2}]]/.Flatten[ans] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 6 :[font = text; inactive; preserveAspect] 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. ;[s] 19:0,0;14,1;21,0;39,1;40,0;44,1;52,0;117,1;121,0;203,1;209,0;288,1;289,0;348,1;349,0;390,1;394,0;424,1;425,0;546,-1; 2:10,13,9,Times,0,12,0,0,0;9,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] 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}]) :[font = text; inactive; preserveAspect] 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. :[font = input; preserveAspect] start=Table[100/3,{i,3},{j,10}]//N; :[font = input; preserveAspect] res=NestList[func,start,500]; :[font = input; preserveAspect] Transpose[res[[500]]]//TableForm :[font = input; preserveAspect; endGroup] Transpose[w]//TableForm :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 7 :[font = text; inactive; preserveAspect] I define the vectors of stay times, mating success and fitness: :[font = input; preserveAspect] 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); :[font = text; inactive; preserveAspect] I define a function to plot fitness against time for a suitable range of values and for a given value of tau : :[font = input; preserveAspect] 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}}]) :[font = text; inactive; preserveAspect] This function can now be used to plot fitness against time for different values of tau, for example: :[font = input; preserveAspect; endGroup] func[4] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 8 :[font = text; inactive; preserveAspect; endGroup] Use pencil and paper. :[font = subsection; inactive; Cclosed; preserveAspect; startGroup] 9 :[font = text; inactive; preserveAspect] I define the male and female frequency distributions and calculate their mean and variance: :[font = input; preserveAspect] 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}; :[font = input; preserveAspect] malemean=Range[31].maledata/Apply[Plus,maledata]//N :[font = input; preserveAspect] femalemean=Range[31].femaledata/Apply[Plus,femaledata]//N :[font = input; preserveAspect] malevar=((Range[31]^2).maledata - maletot malemean^2)/ Apply[Plus,maledata] :[font = input; preserveAspect] femalevar=((Range[31]^2).femaledata - femaletot femalemean^2)/Apply[Plus,femaledata] :[font = text; inactive; preserveAspect] 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]. ;[s] 9:0,0;41,4;46,0;108,4;112,0;195,1;196,0;215,2;218,3;221,-1; 5:4,13,9,Times,0,12,0,0,0;1,14,10,Symbol,0,12,0,0,0;1,18,12,Times,32,10,0,0,0;1,12,9,Times,0,10,0,0,0;2,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] 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])]) :[font = text; inactive; preserveAspect] 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: ;[s] 6:0,0;22,1;34,2;35,0;121,1;129,0;162,-1; 3:3,13,9,Times,0,12,0,0,0;2,13,9,Times,1,12,0,0,0;1,14,10,Symbol,1,12,0,0,0; :[font = input; preserveAspect] Remove[func] func[t_]:=F[t]+f[t]/0.137 :[font = input; preserveAspect] Plot[func[t],{t,0,50}] :[font = input; preserveAspect] tstar=t/.FindRoot[func[t]==1,{t,17.5}] :[font = text; inactive; preserveAspect] 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.) ;[s] 9:0,0;8,1;13,0;48,1;49,0;159,1;164,0;169,1;174,0;210,-1; 2:5,13,9,Times,0,12,0,0,0;4,13,9,Times,1,12,0,0,0; :[font = input; preserveAspect] m[t_]:=If[tt)/0.137,0] :[font = input; preserveAspect] Plot[m[t],{t,0,25}] :[font = input; preserveAspect] mean=NIntegrate[t m[t],{t,0,tstar}] :[font = input; preserveAspect; endGroup; endGroup] variance=NIntegrate[(t-mean)^2 m[t],{t,0,tstar}] ^*)