The cellular automaton model
Cellaut.m
contains matlab code for the cellular automaton model. It is currently set
to run 20 iterations, each of 5000 steps, under the given conditions.
The reaction-diffusion model
This program consists of six
separate .m files. Diffuse.m is the master routine. Diffuse1.m to
Diffuse4.m are diffusion routines that are called by Diffuse.m for
cells with one to four neighbours. Reproduce.m contains the logistic
population increase term. This model provides a framework into which
any set of coordinates can be read. It requires that parameter
restrictions are observed (as outlined in the methods section) and
that the distance between cells is unity. Grids based on real
coordinates can easily be reformatted to give unit coordinate spacing;
simply subtract the minimum x and y coordinates and divide by the
inter-cell distance. If this is done in a spreadsheet, keeping the
key of matching coordinate pairs, the results can readily be matched
to the original coordinates and plotted in ‘real’
space. Note that dispersal can only occur between adjacent cells;
there is no long-distance dispersal to unconnected habitat patches.
In both cases you will need to change the working directory (from
which input files are read, and to which the results file is sent) by
altering the line ‘cd d:\myfolder’ to ‘cd
yourworkingdiskdrive:\yourfolderofchoice’. (e.g. change
‘d:\myfolder’ to d:\models\diffusion’). You may
also need to add your folder to Matlab’s directory using the
‘set path’ option in the Matlab command window.
cellaut.m
% Cellular automaton model - stochastic dispersal, reproduction and mortality
% Written by Graeme Cumming, November 2000
% set the working directory for your own machine
cd e:\cellaut2
% these commands read in x,y coords from file myfile.txt
fid = fopen ('stream0.txt','r');
[xcoord,ycoord] = textread ('stream0.txt','%f %f');
status = fclose (fid);
count = size (xcoord);
% initialise the variables for the program
cell = zeros (count);
cell2 = zeros (count);
cell (500) = 1; %starting value
dispdist = 3; %dispersal distance in m
inc = 2; %number of 'offspring' per iteration
mortality = 300; %death rate / 1000
area = ceil (3.141592654*(dispdist^2)); %max possible number of neighbours
neighbour = zeros (count,area);
reps = 5000;
%first we write a matrix containing identifiers of all cells in dispersal range
for x = 1:count
match = 1;
for i = 1:count
distance = sqrt(((xcoord(x)-xcoord(i))^2)+((ycoord(x)-ycoord(i))^2));
if distance <= dispdist
match = match + 1;
neighbour(x,1) = neighbour(x,1)+1;
neighbour(x,match) = i;
end;
end;
end;
%now begins the main loop of the program
for z = 1:100
cell = zeros (count);
cell2 = zeros (count);
cell (500) = 1;
results = zeros (10000,1);
for r = 1:reps
%disperse
for j = 1:count
if cell(j) == 1
cell2(j) = 1;
for m = 1:inc
rnum = rand * neighbour(j,1);
p = ceil(rnum);
if p==1
p = 2;
end;
plusone_id = neighbour(j,p);
cell2(plusone_id) = 1;
end;
end;
end;
cell = cell2;
cell2 = zeros (count);
for k = 1:count
if round(rand*1000) < mortality
cell(k) = 0;
end;
end;
%results
tally = sum(cell);
results(r,1) = tally;
if tally==0
break
end;
disp (r);
end;
%^ ends the main loop
if z == 1
final = results;
end;
if z > 1
final = cat(2,final,results);
end;
end;
save ('results1.txt','final','-ascii');
diffuse.m
% Program Diffuse
% Written by Graeme Cumming, October 2000
% Program aim is to provide a general reaction-diffusion modelling framework
% Coordinates of study system are read in from a text file 'coords'
% A matrix 'proxim' is then written that holds the number of adjacent cells
% and which cells are adjacent, for each cell
% Diffusion occurs in one of three ways depending on # adjacent cells
% 1 adjacent cell: simple exchange
% 2 adjacent cells: reaction-diffusion in one plane
% 3 adjacent cells: r-d between 3 peripheral and one central cells
% 4 adjacent cells: two dimensional r-d
% After diffusion is completed for each cell, populations increase
% Population increase is determined by a logistic growth equation
% start by setting the home directory
cd e:\diffusn\streams\
% these commands read in x,y coords from columns in file tlake1.txt
fid = fopen ('stream0.txt','r');
[xcoord,ycoord] = textread ('stream0.txt','%f %f');
status = fclose (fid);
count = size (xcoord);
% the coordinate data are now in two vectors, 'xcoord' and 'ycoord'
% 'count' holds the number of coordinates
global adjacent
adjacent = zeros (count,5);
% reserve space for number of and identifiers of adjacent cells
global current;
current = zeros (count);
% the vector 'current' holds current n for each cell at time t
global cplus1;
cplus1 = zeros (count);
% 'cplus1' holds n for each cell at time step t+1
global D;
D = 0.1;
% D is the diffusion coefficient
% for the current diffuse2 it must be <= .5
global cellwidth;
cellwidth = 1;
% the real-world width of the cells, which are presumed to be squares
% now to calculate the proximity key
for j = 1:count
adjacent(j,1) = 0;
for i = 1:count
distance = sqrt(((xcoord(i)-xcoord(j))^2)+((ycoord(i)-ycoord(j))^2));
if distance == cellwidth
adjacent(j,1) = adjacent(j,1) + 1;
celltofill = adjacent(j,1) + 1;
adjacent(j,celltofill) = i;
end;
end;
end;
% the 5 x count matrix has number of neighbours and up to 4 cell identifiers
% for each cell in the analysis
% seed the 500th cell with 100 individuals
current(500) = 100;
% next, the main diffusion loop
iterations = 3000;
step = 1;
%results = current;
num = 1;
results = zeros (count);
for k = 1:iterations
% number of diffusion iterations
% the following calls different diffusion routines
for j = 1:count
switch (adjacent(j,1))
case (1)
diffuse1 (j);
case (2)
diffuse2 (j);
case (3)
diffuse3 (j);
case (4)
diffuse4 (j);
otherwise
disp ('ERROR - not finding neighbours');
end;
reproduce (j);
end;
current = cplus1;
% save the results every n steps
if k/step == floor(k/step)
% results = cat(2,results,cplus1);
results(num) = sum(cplus1);
num = num+1;
end;
disp (k)
if k/3000 == floor(k/3000)
save ('results1.txt','results','-ascii');
end;
end;
diffuse1.m
function [diff1] = diffuse1(j)
% diffusion into the end-cell
global cplus1;
global adjacent;
global current;
global D;
global cellwidth;
neighbour = adjacent(j,2);
unj = current(neighbour);
% unj is the current value of the cell next to the end-point
endcell = current(j);
if endcell <= unj
cplus1(j) = endcell + (unj-endcell)/2;
end;
if endcell > unj
cplus1(j) = endcell - (endcell-unj)/2;
end;
diffuse2.m
function [diff2] = diffuse2(j)
global cplus1;
global current;
global adjacent;
global D;
global cellwidth;
unj = current(j);
currnt1 = adjacent(j,2);
unjplus1 = current (currnt1);
currnt2 = adjacent (j,3);
unjminus1 = current (currnt2);
RHS = (((unjplus1 - (2*unj) + unjminus1))/(cellwidth^2))*D;
% cellwidth = delta x
cplus1(j) = RHS + unj;
diffuse3.m
function [diff3] = diffuse3(j)
global cplus1;
global current;
global adjacent;
global D;
global cellwidth;
unj = current(j);
c1 = adjacent(j,2);
c2 = adjacent(j,3);
c3 = adjacent(j,4);
% c1-c3 are the identifiers of the 3 neighbouring cells
currnt1 = current(c1);
currnt2 = current(c2);
currnt3 = current(c3);
% currnt1-currnt3 are the values of the 3 neighbouring cells
RHS = (((currnt1 - (3*unj) + currnt2 + currnt3))/(cellwidth^2))*D;
% cellwidth = delta x = 1; delta t is 1
% the flow is now three-way
cplus1(j) = RHS + unj;
diffuse4.m
function [diff4] = diffuse4(j)
global cplus1;
global current;
global adjacent;
global D;
global cellwidth;
unj = current(j);
c1 = adjacent(j,2);
c2 = adjacent(j,3);
c3 = adjacent(j,4);
c4 = adjacent(j,5);
% c1-c4 are the identifiers of the 4 neighbouring cells
currnt1 = current(c1);
currnt2 = current(c2);
currnt3 = current(c3);
currnt4 = current(c4);
% currnt1-currnt4 are the values of the 4 neighbouring cells
RHS = (((currnt1 - (4*unj) + currnt2 + currnt3 + currnt4))/(cellwidth^2))*D;
% cellwidth = delta x = 1; delta t is also 1
% the flow is now four-way
cplus1(j) = RHS + unj;
reproduce.m
function [shithappens] = reproduce (j)
% exponential population increase
global cplus1;
N = cplus1(j);
r = .5;
K = 500;
babies = r*N*(1-(N/K));
%fate = rand;
% fate is a random number between 0 and 1
% corpses = ((fate/1)*N)/2;
% a random proportion of the population, not exceeding 50%, dies at every iteration
corpses = 0;
cplus1(j) = cplus1(j) + babies - corpses;