function D = fsim ( Npops, mu0, mu1, gamma0, gamma1, IC, NReps, MaxTime )
%function to perform metapopulation simulations
%Graham Medley, February 2013; University of Warwick
%data structure for the simulation output. Initially this is a 3-D array
%consisting of output at the start of each time interval (and the end of
%the last one) for each repetition and each population
D = zeros ( NReps, MaxTime+1, Npops );
for irep = 1 : NReps
pops = IC; %pops is set to the initial condition
time = 0; %time is set to zero
outcount = 1; %this is a counter to record the time interval
D ( irep, outcount, : ) = pops; %record the initial condition as the first output
outtime = 1; %this is the time of the next output
%perform a translocation ( at the start )
pops = translocate ( time, pops );
%continue until MaxTime or metapopulation extinction or one population
%gets silly
while ( time < MaxTime ) && any ( pops > 0 ) && all ( pops < 100 )
BirthRates = max ( 0, gamma0 - gamma1 .* pops );
DeathRates = mu0 + mu1 .* pops;
BR = pops .* BirthRates; %calculate the population level birth rates
DR = pops .* DeathRates; %calculate the population level death rates
EventRates = cumsum ( [ sum(BR) sum(DR) ] );
TotalRate = EventRates ( end ); %the total rate of events
randnums = rand ( 1, 3 ); %get 3 random numbers from a uniform distribution btn. 0 and 1
tstep = - log( randnums(1) ) / TotalRate; %calculate the time to the next event using the 1st random number
time = time + tstep; %time moves ever on...
oldpops = pops; %remember the population in its current state
%this section chooses the event and the population, and changes the
%population accordingly
%use the second random number to decide btn. events
eventindex = find ( EventRates > randnums ( 2 ) * TotalRate );
switch eventindex ( 1 )
case 1
PopRate = cumsum ( BR );
index = find ( PopRate > ( randnums ( 3 )*PopRate ( end ) ) );
pops ( index ( 1 ) ) = pops ( index ( 1 ) ) + 1; %add one to the population
case 2
PopRate = cumsum ( DR );
index = find ( PopRate > ( randnums ( 3 )*PopRate ( end ) ) );
pops ( index ( 1 ) ) = pops ( index ( 1 ) ) - 1; %remove one from the population
otherwise
%this indicates an error - report it and stop simulating
disp ( 'error in fsim' );
return
end
while ( time > outtime ) && ( outtime <= MaxTime )
outcount = outcount + 1; %record the next time point
D ( irep, outcount, : ) = oldpops; %population as it is at outtime
outtime = outtime + 1; %next time point for recording
pops = translocate ( time, pops ); %do the translocation annually
end
end %while
end %for NReps