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

In summary, the conversation revolved around a code for converting biomass of elk into new wolves using ODE and Gillespie methods. The code includes a reversible reaction and simulations for both methods. The speaker also mentioned plans to upload several hundred simple mfiles for public use.
  • #1
DrWahoo
53
0
Copy and paste if you want to use any of this code. Its a simple iteration of converting biomass of elk into new wolves.

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');
View attachment 7493
 

Attachments

  • MHB.png
    MHB.png
    3 KB · Views: 71
Physics news on Phys.org
  • #2
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.
 
  • #3
MarkFL said:
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.
 
  • #4
DrWahoo said:
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)
 

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

1. How does the Gillespie algorithm work in the context of a reaction model with wolves and elk?

The Gillespie algorithm is a computational method used to simulate stochastic processes in systems with discrete events, such as chemical reactions. In the context of a reaction model with wolves and elk, the algorithm can be applied to simulate the population dynamics of both species based on their interactions and birth/death rates.

2. What are the key parameters that need to be included in the Matlab code for this reaction model?

The key parameters that need to be included in the Matlab code for this reaction model are the birth/death rates of wolves and elk, the carrying capacity of the environment, and the interaction rates between the two species. Additionally, the initial population sizes of wolves and elk and the simulation time step should also be specified.

3. How can I validate the accuracy of the simulation results from the Matlab code?

One way to validate the accuracy of the simulation results is to compare them with real-world data or with results from other established simulation models. Sensitivity analysis can also be performed by varying the input parameters and observing the changes in the simulation output.

4. Are there any limitations or assumptions in using the Gillespie algorithm for this reaction model?

One potential limitation of using the Gillespie algorithm for this reaction model is that it assumes well-mixed populations and does not account for spatial dynamics. Additionally, it may not accurately capture rare events or large fluctuations in population size due to its stochastic nature.

5. Can the Matlab code be modified to include other factors or species in the reaction model?

Yes, the Matlab code can be modified to include other factors or species in the reaction model. The key parameters and equations would need to be adjusted accordingly to account for the new additions. However, it is important to ensure that the model remains biologically realistic and that the simulation results are still valid.

Similar threads

  • MATLAB, Maple, Mathematica, LaTeX
Replies
6
Views
2K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
2
Views
1K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
4K
  • Engineering and Comp Sci Homework Help
Replies
4
Views
3K
Back
Top