Matlab Sim Code for ODE & Gillespie for Reversible Reaction

[DrWahoo]

Here is the code you input into matlab. Aini etc, are the initial values of the population densities. A for predator, B for Prey.

% 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');

 

[Joppy]

Bit hard to read...

 

[DrWahoo]

I am trying to upload a PDF of what the chemical reaction is doing to help the logic in the iterations. Maybe Mark can help me upload it. Its saying my file size is to large.

 

[Joppy]

This was part of an exercise I did sometime last year I think. It is larger than it needs to be since I was solving the equations with a manually implemented RK4 scheme.

EDIT: Wow. The forum recognizes abbreviations, and I didn't even have to do anything!

 

[S.G. Janssens]

I think the code contributions in this and other recent threads are useful and unselfish, which is nice.

Would it be an idea to set up a small MHB code repository (perhaps using git or svn)? This would allow for easier forking, editing and keeping track of authors and changes. (In short, it would help to maintain an overview.)

Of course, I cannot oversee how feasible this is, but the possibility crossed my mind.

 

[MarkFL]

I am trying to upload a PDF of what the chemical reaction is doing to help the logic in the iterations. Maybe Mark can help me upload it. Its saying my file size is to large.

Did you try zipping the file before uploading it?

 

[DrWahoo]

Did you try zipping the file before uploading it?

Yeah, still won't upload. Its 3.1 Megs.
I will shoot you an email Mark.