% Lotka model - uses deterministic and stochastic ODEs
function results = LotkaStoch
% Wendy Thomas 2010
% solves Lotka Voltare predator prey model:
% H --> 2H, rate r (prey reproduction)
% H + P --> 2*b/a*P, rate a (predator reproduction and prey death)
% P --> nothing, rate m (predator death)
clear all; close all;
r = 2;
a = 0.01;
b = 0.01;
m = 1;
H0 = 100;
P0 = 200;

IC = [H0 P0];
tspan = [0,20];
params = [r,a,b,m];
figure(1)
% deterministic model:
[T,Y] = ode23s(@LotkaODE,tspan,IC,[],params);
plot(T,Y); hold on;
results.det = [T,Y];
% stochastic models:
[t,y] = DSDEexact(@LotkaRXN,tspan,IC, params);
plot(t,y,'--');
results.esa = [t,y];
[ttau,ytau] = DSDEtauleap(@LotkaRXN,tspan,IC, 0.01, params);   
plot(ttau,ytau, ':');
results.tau = [ttau,ytau];


function dydt = LotkaODE(t,y, params)
r = params(1);
a = params(2);
b = params(3);
m = params(4);
H = y(1); P = y(2);  
dH = r*H-a*H*P;
dP = b*H*P-m*P;
dydt = [dH; dP];


function [S, P, K] = LotkaRXN(params)
%defines the reactions and the stochastic rate constants for:
% H --> 2H, rate r
% H + P --> 2*b/a*P, rate a
% P --> nothing, rate m
r = params(1);
a = params(2);
b = params(3);
m = params(4);
% Need a row for each reaction and a column for each chemical species
%   H P
S = [  % How many of each chemical is used as substractes
    1 0
    1 1
    0 1]  ;
P = [  % How many does each chemical species is made as products
    2 0
    0 2*b/a 
    0 0] ;
K  = [r
    a
    m];


