% SIR model - uses deterministic and stochastic ODEs
function results = SIR
% Wendy Thomas 2010
% solves susceptible-infectious-resistant disease model:
%  S + I --> 2I with rate constant beta
%      I --> R with rate constant nu
clear all; close all;
beta = 0.001;
nu=0.3;
S0 = 1000;
I0= 2;
tspan = [0,40];
ratio = beta/nu*S0 % if > 1, outbreak
IC = [S0,I0,0];
params = [beta,nu];
% deterministic model:
[T,Y] = ode23s(@SIRode,tspan,IC,[],params);
plot(T,Y); hold on;
results.det = [T,Y];
% stochastic models:
[t,y] = DSDEexact(@SIRrxns,tspan,IC, params);
plot(t,y, ':');
results.esa = [t,y];
[ttau,ytau] = DSDEtauleap(@SIRrxns,tspan,IC, 0.1, params);   
plot(ttau,ytau,'-.');
results.tau = [ttau,ytau];


function dydt = SIRode(t,y, params)
beta = params(1); nu = params(2);
S = y(1); I = y(2);  R = y(3);
dS = -beta*I*S;
dI = beta*I*S - nu*I;
dR = nu*I;
dydt = [dS;dI;dR];


function [S, P, K] = SIRrxns(params)
%defines the reactions and the stochastic rate constants for:
%  S + I --> 2I beta
%      I --> R nu
% Need a row for each reaction and a column for each chemical species
%   S  I  R
S = [  % How many of each chemical is used as substractes
    1  1  0
    0  1  0]  ;
P = [  % How many does each chemical species is made as products
    0  2  0
    0  0  1] ;
K  = params';
% parameters will be passed separately as beta, nu


