% Written by E de los Santos for BE150, Winter 2013

% This script implements a simple euler integration for use with delayed
% differential equations. The script requires the following inputs:
 
% model - a function that takes in 4 arguments (current mRNA level,
% current protein level, previous mRNA level, previous protein level). It
% should return a two dimensional output [dm/dt ; dp/dt]

% delaym - the length of the delay of the mRNA level (Tp)

% delayp - the length of the delay of the protein level (Tm)

% tstep - size of the time step, should be small enough

% tend - end time for the integration (we are assuming we start at t=0)

% It returns a time vector and the protein and mRNA levels for each time
% step in the time vector.


function [t,m,p] = dde(model, delaym, delayp, tstep,tend),
t(1) = 0;
p(1) = 0;
m(1) = 0;

ndelaym = floor(delaym/tstep);
ndelayp = floor(delayp/tstep);
nfinal = floor(tend/tstep);  

for i=2:nfinal,
    
    if (i > ndelaym),
        mdelay = m(i-ndelaym);
    else
        mdelay = 0;
    end
    
    if (i > ndelayp),
        pdelay = p(i-ndelayp);
    else
        pdelay = 0;
    end
    
    dx = model(m(i-1),p(i-1),mdelay,pdelay);
    m(i) = m(i-1) + tstep*dx(1);
    p(i) = p(i-1) + tstep*dx(2);
    
    t(i) = i*tstep;
end

end
