
MBW:Agent Based Model of Fish PopulationsFrom MathBioContentsExecutive SummaryConstructing ecological models can be challenging. There are many complex interactions between species, and building a model to understand the behavior of the ecological system is not usually trivial. Agent based models are one tool to understand ecological and other complex systems. Agent based models can easily account for heterogeneity in a population. The technique is a Monte Carlo simulation that tracks individuals in a population in discrete time. James Rice and Thomas Miller in their 1993 paper, “Growth Rate Variation and Larval Survival: Inferences from an IndividualBased SizeDependent Predation Model” construct an agent based model to understand larval bloater populations in response to predation from alewife. They explore how mean and variance of the initial larval bloater population responds to predation. Similar agent based models have proved useful in a management context, particularly when crossscale dynamics (individuals to populations) are complex (e.g. inSTREAM). Context/Biological situationThe biological system is that of a body of water, where the alewife fish (Alosa pseudoharengus) preys on larval/juvenile bloater (Coregonus hoyi). An example habitat for the two species is Lake Michigan. The model is designed to investivage the connection between growth rates of bloater and selective pressure on the bloater population. Some studies suggest that small decreases in initial growth rates cause a large increase in total mortality, because the population will be most vulnerable longer. Spread in diversity of sizes promotes survivability, and mean growth rates will be higher over time. The simulation takes place over a period of 60 days. The model assumes that the growth rate of bloater is much higher than growth rate of alewife. This allows us to hold alewife size constant in the simulation. This is reasonable, as alewife start out more mature, and grow slowly. Two hypotheses about the biological system are being tested. The first hypothesis is that increasing either the mean growth rate or the variability in growth rate will lead to an increase in survival rate when predation is based on size, and that the survivors will tend to have higher growth rates than average. The second hypothesis tests the assertion that a population with growth rates that change from day to day in a random walk fashion will have different size distribution and survival rate than a population where the growth rates of individual agents are held constant though the course of the simulation. ModelThe model is agent based: Bloater fish will be tracked individually by the software, over a sixty day period. Alewife are not treated individually. Their interactions with bloaters are statistical, based on previous field studies. This enters the model as a product of three probabilities: chance of encountering an alewife, chance of being attacked, and chance of being eaten if attacked. Probability of being eaten if attacked is the only one of these parameters that has been measured. The others are educated guesses. Encounter rate is held constant, because this is a function of swim speed of the predator. The speed of the bloater is much less than alewife speed, so it is not very important. It is also known the probability of being eaten if attacked drops sharply with bloater size. The chosen value of contact probability was 0.2. This is less than the actual value, however the actual values causes the population to collapse before the end of the 60 days. A value of 0.2 allows the simulation to run, but still have a high enough mortality rate to investigate changes in population. The model also assumes there will not be multiple encounters in one days. The bloaters start out all the same size, 12mm and the alewife are all 90mm. The growth rate per day is set at the beginning, according to a normal distribution. In the simulation we vary the mean between 0.2, 0.4 and 0.6 mm/day. The standard deviation also varies between 0.04, 0.08 and 0.16 mm/day. We do one simulation where the growth rate changes randomly, with the mean moving as a random walk. For the nonrandom growth rates, we clamp the maximum rate to 1mm per day, and minimum at 0 mm/day. For the random walk growth rates, we clamp the maximum rate to .8mm per day, and minimum at 0 mm/day to make it symmetrical around the starting growth rate of .4mm/day. To initiate a day, we first generate a random number between 0 and 1. If this number is less than contact probability, we run an equation to determine if the fish will die or not. This equation is:
We then generate a second random number. If this number is less than the probability of dying, the fish will die. Built into this equation is if bloater is longer than 40.9mm, it will not be eaten. After the deaths are calculated, the survivors are grown according to their initial growth rates. Simulations were run with between 4000 and 10000 individual bloaters. The MATLAB code for the agent can be found in the Appendix. For more explanation, and to see another application, of predatorprey models, see MBW:PredatorPrey Models With Delay And Prey Harvesting. AnalysisWe ran multiple ensembles of simulations to understand the behavior of populations. We wanted to understand how mean growth rate and standard deviation in growth rate affected the population at the end of 60 days. The authors made frequency distributions showing how these variables changed the initial and final populations.
The red is the final, and blue initial. The x axis is the size of the bloater. The y axis is the percent of the current population. From the data, it can be seen that mean growth rate has a profound effect on survivability, and variance has an effect, but it is not as strong. A higher variance shifts the average mean growth rate to the right. This is higher growth rates being selected for over lower growth rates. The largest change in distribution is a low mean growth rate with high variance. In this case mean growth rate is not overpowering the effects of variability. We also look at what percent of the initial cohort survived. Mean growth rate has the strongest impact on these values. Below is a table of final mean growth rates, ordered by initial mean growth rate and standard deviation. The order is the same as the 3x3 grid in the previous figure.
In addition to growth rates we looked at a frequency distribution of final lengths (in mm). These show a similar pattern as the growth rate distributions, with high variability causing a shift to the right. This indicates that survivors are selected for based on size. Now that the initial and final distributions are known, how do the populations behave over the 60 day period? We ran the model and looked at variables as functions of time, or variables parameterized by time, so we can understand the dynamics. First is the number of survivors as a function of time. It is apparent that this is a nonlinear relationship. The bloaters die rapidly at first, and their population stabilizes near the end of 60 days. This is indicating they are dying more slowly, because they are growing larger. This can be seen in a different way, by looking at mean growth rate as a function of total mortality. The mean growth rate is approximately 0.2 mm/day until 60% of the population has died. At this point the mean growth rate rapidly increases. This is indicating that predation is uniform in the population until most of the cohort is dead, at which point mean growth rate increases rapidly, favoring the stronger bloaters. We wanted to see what growth rate as a function of time, which the authors did not publish. This figure indicates that growth rates are indeed initially slow. All of the simulations are flat at small time, and have a large approximately linear region. The authors made some theoretical predictions about the probability of surviving to 60 days, given how long they have already survived. They also calculated the daily probability of dying given survival to some time. They did not compare their simulation with the theory, however we did. This graph was generated with the probability of dying equation above. The two above equations were generated in MATLAB according to the formula The experiment fits closely on the theory line, except at low growth rates. For higher standard deviations, the theory and experiment begin to strongly disagree. This may indicate a flaw in the model. In a second simulation the authors looked at the effect of a nonconstant growth rate. At each day, a new growth rate for an individual is selected. The new growth rate is selected from a normal distribution of SD = 0.8, and a mean of the present growth rate. If the generated growth rate is above 0.8, or less than 0, it is set to 0.8 or 0. We ran this model with and without predation. Without predation the variability increased vs constant growth, as expected. With predation, like in previous simulations, there was a shift of the mean to higher growth rates. The cohort survival also increased 2.81% to 3.23%. We also decided to look at the effect of predator contact probability on surviving population (out of 10,000), Mean growth rate vs predator contact, and finally surviving population vs mean growth rate: ConclusionsThe results indicated that changes in variance of growth rate and mean growth rate do have effects on the differences in mean growth rate and mean size between the initial population and the final population of survivors. Changes in these two parameters affect the survival percentage as well. Our results indicated that increasing the mean growth rate and variance will lead to increases in mean size and mean growth rate of survivors. Survival rates were also increased by increasing mean growth rate, and in most but not all cases increasing growth rate variance had a similar effect. These results support all of the assertions in the first hypothesis except that increasing growth rate variance will lead to an increase in survival. An increase in variance was associated with an increase in survival when mean growth rate was low, but with higher mean growth rate increases in growth rate variance didn't always show an increase in survivalas was seen when mean growth rate was equal to .4 and variance was increased from 0.04 to 0.08. The results also supported the second hypothesis. With mean growth rate equal to 0.4 and standard deviation of growth rate equal to 0.08, simulating the random walk model 30 times yielded a mean survival percentage of 3.23, whereas 30 simulations of the non random walk model yielded a mean survival percentage of 2.81. This supports the assertion that a population with random walk effects in the mean growth rate experience higher survival rates. The data also showed that the random walk effect increased the variation in size distribution of the survivors when compared to the model without the random walk effects. AppendixWe ran all the simulations in Matlab R2010a.
Here's the code for the random walk ABM: clear g x = [0:.05:1]; nfish=100000; ndays=60; %100000 fish for 60 days sizes=zeros(ndays,nfish); growth=zeros(ndays,nfish); alive=ones(1,nfish); sigmag=0.08; % standard deviation of changes in growth % initial sizes and growth rates sizes(1,1:nfish)=12*ones(1,nfish); growth(1,1:nfish)=0.4; %iterate the model for jday=1:(ndays1); % Predation. For simplicity of coding, the dead can die again. meetpred=rand(1,nfish)<0.7; %does a predator find me? pdie=0.33 + 0.15*(90./sizes(jday,:)); %if so.... eaten=rand(1,nfish)<pdie; m=find(meetpred.*eaten>0); % do I live or die? alive(m)=0; % find sizes and growth rates for tomorrow sizes(jday+1,:)=sizes(jday,:)+growth(jday,:); growth(jday+1,:)=growth(jday,:)+sigmag*randn(1,nfish); growth(find(growth<0))=0; growth(find(growth>0.8))=0.8; m=find(alive>0); growthrate=(sizes(jday+1,m)sizes(jday,m)); g(jday,:) = hist(growthrate,x)./sum(m); end; The code to create figure 1 is two files, an agent and a script to run it. Script: gr = [0.2, 0.4, 0.6]; sd = [0.04, 0.08, 0.16]; k=1; for y=1:3 for x=1:3 subplot(3,3,k); agent(gr(x), sd(y)) table1(y,x,:) = runragent(gr(x), sd(y)); k = k + 1; end end Agent: function ps = agent(ms, ss) x = [0:.05:1]; nfish=4000; ndays=60; %5000 fish for 50 days sizes=zeros(ndays,nfish); growth=zeros(ndays,nfish); alive=ones(1,nfish); %sigmag=sig; % standard deviation of changes in growth % initial sizes and growth rates sizes(1,1:nfish)=12*ones(1,nfish); growth(1,:)= ms + ss*randn(1,nfish); %G; growth(find(growth(1,:)<0))=0; growth(find(growth(1,:)>1.0))=1.0; % h=hist(growth(1,:), x)/nfish; % bar1=bar(x,h,'Facecolor','b','EdgeColor','b'); % set(bar1,'BarWidth', 1.0); % title('Initial'); %iterate the model for jday=1:(ndays1); % Predation. For simplicity of coding, the dead can die again. meetpred=rand(1,nfish)<0.2; %does a predator find me? pdie=0.33 + 0.15*(90./sizes(jday,:)); %if so.... eaten=rand(1,nfish)<pdie; m=find(meetpred.*eaten>0); % do I live or die? alive(m)=0; % find sizes and growth rates for tomorrow sizes(jday+1,:)=sizes(jday,:)+growth(1,:);%growth(jday,:); %growth(jday+1,:)=growth(jday,:) + sigmag*randn(1,nfish); %random walk %growth(find(growth<0))=0; %growth(find(growth>0.8))=0.8; m=find(alive>0); growthrate=(sizes(jday+1,m)sizes(jday,m)); g(jday,:) = hist(growthrate,x)./sum(m); end; % plot growth rate of survivors % figure; m=find(alive>0); growthrate=(sizes(ndays,m)sizes(1,m))/ndays; % hold on; h=hist(growthrate,x)/length(m); bar2=bar(x,h,'FaceColor','b','EdgeColor','b'); set(bar2, 'BarWidth',0.7); % hold off; % title(num2str(length(m)/nfish*100)); %figure(2); %surf(x(1,:),1:jday,g); %lv=mean(growthrate); %rv=std(growthrate); ps = length(m)/nfish*100; Code for generating theory graphs: d=60; % Pd = 066 + 0.03*(90/12) grl = [0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0]; for g=1:7 gr = grl(g); for t=0:d Pd(g,t+1) =( 0.066 + 0.03*(90/(12 +gr*t)) ); if(Pd(g,t+1) < 0) Pd(g,t+1) = 0; end end end % plot(tx,Pd(1,:), tx, Pd(2,:), tx, Pd(3,:),tx, Pd(4,:),tx, Pd(5,:),tx, Pd(6,:),tx, Pd(7,:) ) % Pd( find(Pd < 0) ) = 0; for c=2:6 for k=1:d for t=k:d if (t==k) Ps(c1,k) = 1Pd(c,t); else Ps(c1,k) = Ps(c1,k)*(1Pd(c,t)); end end end end tx = 1:length(Ps(1,:)); plot(tx, Ps(1,:),tx, Ps(2,:),tx, Ps(3,:),tx, Ps(4,:),tx, Ps(5,:) ) hold on; SD=0.16; MN = 100; for c=1:5 gr = grl(c+1); Pt = agentPgs(gr, SD); for e=2:MN Pt = Pt + agentPgs(gr, SD); end Psd(c,:) = Pt/MN; plot(tx, Psd(c,:),'o'); end hold off; Associated Agent: function ps = agentPgs(ms, ss) x = [0:.05:1]; nfish=4000; ndays=60; %5000 fish for 50 days sizes=zeros(ndays,nfish); growth=zeros(ndays,nfish); alive=ones(1,nfish); %sigmag=sig; % standard deviation of changes in growth % initial sizes and growth rates sizes(1,1:nfish)=12*ones(1,nfish); growth(1,:)= ms + ss*randn(1,nfish); %G; growth(find(growth(1,:)<0))=0; growth(find(growth(1,:)>1.0))=1.0; % h=hist(growth(1,:), x)/nfish; % bar1=bar(x,h,'Facecolor','b','EdgeColor','b'); % set(bar1,'BarWidth', 1.0); % title('Initial'); %iterate the model for jday=1:(ndays1); % Predation. For simplicity of coding, the dead can die again. meetpred=rand(1,nfish)<0.2; %does a predator find me? pdie=0.33 + 0.15*(90./sizes(jday,:)); %if so.... eaten=rand(1,nfish)<pdie; m=find(meetpred.*eaten>0); % do I live or die? alive(m)=0; % find sizes and growth rates for tomorrow sizes(jday+1,:)=sizes(jday,:)+growth(1,:);%growth(jday,:); %growth(jday+1,:)=growth(jday,:) + sigmag*randn(1,nfish); %random walk %growth(find(growth<0))=0; %growth(find(growth>0.8))=0.8; m=find(alive>0); growthrate=(sizes(jday+1,m)sizes(jday,m)); if jday==1 ns(jday,:) = m; else ns(jday,:) = [m,zeros(1, length(ns(1,:))length(m) )]; end % g(jday,:) = hist(growthrate,x)./sum(m); end; % plot growth rate of survivors % figure; m=find(alive>0); growthrate=(sizes(ndays,m)sizes(1,m))/ndays; % hold on; % h=hist(growthrate,x)/length(m); % bar2=bar(x,h,'FaceColor','r','EdgeColor','r'); % set(bar2, 'BarWidth',0.7); % hold off; % title(num2str(length(m)/nfish*100)); %figure(2); %surf(x(1,:),1:jday,g); %lv=mean(growthrate); %rv=std(growthrate); for j=1:(ndays1) s1=find( ns(j,:) > 0); s2 = find( m > 0); a=intersect( s1,s2); if length(s1) == 0 fr(j) = 0; else fr(j) = length(a)/length(s1); end end fr(ndays)=1.0; ps = fr; Code for generating functions of predation rate: for k = 1:30; p = changingpredation(k/43); r(k) = p(1); y(k) = p(2); end; And the agent: function rv = changingpredation(predrate) x = [0:.05:1]; nfish=10000; ndays=60; %10000 fish for 60 days sizes=zeros(ndays,nfish); growth=zeros(ndays,nfish); alive=ones(1,nfish); sigmag=0.08; % standard deviation of changes in growth % initial sizes and growth rates sizes(1,1:nfish)=12*ones(1,nfish); growth(1,1:nfish)=0.4; %iterate the model for jday=1:(ndays1); % Predation. For simplicity of coding, the dead can die again. meetpred=rand(1,nfish)< predrate; %does a predator find me? pdie=0.33 + 0.15*(90./sizes(jday,:)); %if so.... eaten=rand(1,nfish)<pdie; m=find(meetpred.*eaten>0); % do I live or die? alive(m)=0; % find sizes and growth rates for tomorrow sizes(jday+1,:)=sizes(jday,:)+growth(jday,:); growth(jday+1,:)=growth(jday,:)+sigmag*randn(1,nfish); growth(find(growth<0))=0; growth(find(growth>0.8))=0.8; m=find(alive>0); growthrate=(sizes(jday+1,m)sizes(jday,m)); g(jday,:) = hist(growthrate,x)./sum(m); end; rv = [length(m), mean(sizes(m))]; Paper SummaryMathematics Used: Monte Carlo simulation Type of Model: individual/agent based Biological System Studied: Dynamics of fish populations based on the behavior of individual fish.
Subsequent ResearchAgent based models such as the one presented here have gained much traction in fisheries research. The utility of an agent based approach rather than a populationbased approach, is that individual behavior can be accounted for explicitly, and it is relatively easy to incorporate variation in individual behaviors. In part because of this feature, agent based modeling has undergone a rapid diversification in the fields of ecology and evolution, as reviewed by DeAngelis and Mooij (2005). DeAngelis and Mooij (2005) detail applications of agent based models in conservation biology, evolution, metapopulation theory, and forest management  all fields where scaling up from individuals (whether they be genes, animals, or trees) to higher levels of organization (i.e. populations, forests) are critical. For an example in community ecology, see MBW:Habitat heterogentiyarea tradeoffs and species diversity. 