- #1

DrWahoo

- 53

- 0

Code:

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