# Matlab Sim Code for ODE & Gillespie for Reversible Reaction

• MATLAB
• DrWahoo
In summary, the population densities of predators and prey are initialized to Aini and Bini, respectively. A reversible reaction takes place between the two populations, with b1 = .033 and bo = .00047. The resulting population densities are A for predator, B for Prey.
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.
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');

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.

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!

#### Attachments

• PredatorPreyExample.zip
2.5 KB · Views: 67
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.

DrWahoo said:
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.

MarkFL said:

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

## 1. What is Matlab Sim Code for ODE for Reversible Reaction?

Matlab Sim Code for ODE (Ordinary Differential Equation) for Reversible Reaction is a method used to simulate the behavior and dynamics of a reversible chemical reaction using Matlab software. It involves solving a set of differential equations that describe the rates of change of the reactants and products in the reaction.

## 2. How does Matlab Sim Code for ODE for Reversible Reaction work?

The Matlab Sim Code for ODE for Reversible Reaction works by numerically solving the set of differential equations using the built-in functions and solvers in Matlab. These equations describe the changes in the concentrations of the reactants and products over time, taking into account the rate constants and stoichiometry of the reaction.

## 3. What is Gillespie Algorithm for Reversible Reaction?

The Gillespie Algorithm for Reversible Reaction is a stochastic simulation method used to model and analyze the kinetics of a reversible chemical reaction. It is based on the concept of randomly selecting and executing individual reaction events, taking into account the reaction rates and probabilities at each step.

## 4. How is Gillespie Algorithm for Reversible Reaction implemented in Matlab?

In Matlab, Gillespie Algorithm for Reversible Reaction can be implemented by writing a custom code using the built-in functions and random number generators. Alternatively, there are also open-source Matlab toolboxes available that have pre-written functions for implementing the Gillespie Algorithm for Reversible Reaction.

## 5. What are the advantages of using Matlab Sim Code and Gillespie Algorithm for Reversible Reaction?

The use of Matlab Sim Code and Gillespie Algorithm for Reversible Reaction provides a more accurate and detailed understanding of the dynamics and behavior of a reversible reaction compared to traditional analytical methods. It also allows for the simulation of complex reaction networks and the incorporation of noise and variability in the reaction system.

• MATLAB, Maple, Mathematica, LaTeX
Replies
3
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
5K
• MATLAB, Maple, Mathematica, LaTeX
Replies
2
Views
1K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
7
Views
3K
• MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
4K
• MATLAB, Maple, Mathematica, LaTeX
Replies
4
Views
2K
• MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
7K
• MATLAB, Maple, Mathematica, LaTeX
Replies
7
Views
3K