MATLAB Matlab code for a Gillespie type reaction (model) with wolves and elk.

[DrWahoo]

Copy and paste if you want to use any of this code. Its a simple iteration of converting biomass of elk into new wolves.

% example for ODE and Gillespie

% one reversible reaction
 b1 = .033;
bo = .00047;
a1 = .0022;
ao = .00055;

Aini = 5168;
Bini = 34;

%% Basic ODE simulation
dt = 0.01;
Tmax = 1;

TV = (0:dt:Tmax)';
NT = length(TV);

ABVals = zeros(size(TV,1),2);

A=Aini;
B=Bini;
ABVals(1,:) = [A, B];

% loop over times
for it=2:NT
    dAdt = - b1*A*B + bo*A;
    dBdt = a1*B*A -ao*B 
    A = A + dAdt * dt;
    B = B + dBdt* dt;
    
    ABVals(it,:) = [A, B];
    
    
end

figure(1)
clf;
plot(TV,ABVals);
xlabel 'T';
legend('A','B');
hold on

return;
%% Gillespie

% loop goes over time steps

%Tmax = 1;
MaxSteps = 1000;

% define a large array to record everything
ABRecord = zeros(MaxSteps,2);
TRecord = zeros(MaxSteps,1);ABCount = [Aini,Bini];

t=0;TRecord(1) = t;
ABRecord(1,:) = ABCount;
istep=2;
while (t<= Tmax  && istep<=MaxSteps),
    
    A=ABCount(1);
    B=ABCount(2);
    
    % propensities ('gamma')
    PropElkEaten = b1 * A*B;
    PropWolvesBorn = a1*A*B;
    PropWolvesDie = ao*B;
    PropElkBorn = bo*A;
    
    % distribution of next reaction times follows exp( - gamma * t)
    NextFwdTime = exprnd(1/PropElkEaten,1);
    NextRevTime = exprnd(1/PropWolvesBorn,1);
    NextWolvesDie = exprnd(1/PropWolvesDie,1);
    NextElkBorn = exprnd(1/PropElkBorn,1);
    
    if NextFwdTime < NextRevTime,
        % forward reaction (A->B) happens
        t = t + NextFwdTime;
        ABCount = ABCount + [-1,1];
    else
        % reverse reaction (A->B) happens
        t = t + NextRevTime;
        ABCount = ABCount + [1,-1];
    end
    
    % record
    ABRecord(istep,:) = ABCount;
    TRecord(istep) = t;
    
    % done with this step, increment the step count
    istep=istep+1;
    
end

if istep-1<MaxSteps,
    TRecord = TRecord(1:istep-1);
    ABRecord = ABRecord(1:istep-1,:);
end

figure(1)
stairs(TRecord,ABRecord)
legend('A - ODE','B - ODE','A - Gillespie','B - Gillespie');

View attachment 7493

 

[MarkFL]

I've enclosed your code in our [CODE][/CODE] tags to preserve the white space, and to make it easier for people to select and copy the entire block by double-clicking it.

 

[DrWahoo]

I've enclosed your code in our [CODE][/CODE] tags to preserve the white space, and to make it easier for people to select and copy the entire block by double-clicking it.

Thanks Mark!

I have several hundred simple mfile's I plan to upload for people to use.

 

[MarkFL]

Thanks Mark!

I have several hundred simple mfile's I plan to upload for people to use.

That's awesome! That's definitely music to my ears! (Yes)