function RiceiGEMmodel2017()
%Rice iGEM Main Model
%Jack Terrell
%1 November 2017
pb = gen_bioc_params()
ps = gen_system_params()
psim = gen_sim_params('tspan',[0,250],'purges',[[48:24:240]',ones(9,1)*120000],'seeds',[[48:24:240]'+0.1,ones(9,1)*150])
[tout,yout] = simulate_system(pb,ps,psim);
web('https://www.youtube.com/watch?v=ZcDj46elLhQ')
end

function [tout,yout] = simulate_system(pb,ps,psim)
%vector specification
%   Aggregate Engineered Cell Volume (L)
%   Aggregate Non-Engineered Cell Volume (L)
%   BOD (mg O_2/L)
%   Cr6_ext (uM)
%   Cr6_int_eng (uM)
%   Cr6_int_nat (uM)
%   LacI (uM)
%   BamHI (uM)
%   chr6 (uM)

%   flowfun (L/hr)
%   bod_in_fun (mg O2/L)
%   purges (purge table) (nx2 matrix) (times of purges (hr), volumes of purges (L))

ddt = @(t,v)(system_derivative(t,v,pb,ps,psim));
tcur = psim.tspan(1);
pidx = 1;
sidx = 1;
tout = [];
yout = [];
v0 = zeros(1,9);
v0(1) = 0.03;
v0(2) = 0.10;
yout(1,:) = v0';
tout(1,:) = psim.tspan(1);
opts = odeset('AbsTol',1e-6,'RelTol',1e-6);
while true
    if sidx>size(psim.seeds,1)&&pidx>size(psim.purges,1)
        tint = psim.tspan(2);
        [tt,yt] = ode15s(ddt,[tcur,tint],yout(end,:)',opts);
        yout = [yout;yt(2:end,:)];
        tout = [tout;tt(2:end,:)];
        break
    elseif sidx>size(psim.seeds,1)
        tint = psim.purges(pidx,1);
        if tint~=tcur
            [tt,yt] = ode15s(ddt,[tcur,tint],yout(end,:)',opts);
            yend = yt(end,:)';
            tend = tt(end);
            yt([end],:) = interp1(tt,yt,[tt(end)-eps(tt(end))]);
            tt(end) = tt(end)-eps(tt(end));
            yout = [yout;yt];
            tout = [tout;tt];
        else
            yend = yout(:,end)';
            tend = tout(end);
        end
        y = yend;
        y(1) = y(1)-(psim.purges(pidx,2))*(y(1)/(y(1)+y(2)))*((1/ps.SVI)/(1-ps.watfrac));
        y(2) = y(2)-(psim.purges(pidx,2))*(y(2)/(y(1)+y(2)))*((1/ps.SVI)/(1-ps.watfrac));
        if tint~=tcur
            yout = [yout;y'];
            tout = [tout;tint];
        else
            yout(end) = y';
        end
        tcur = tint;
        pidx = pidx+1;
    elseif pidx>size(psim.purges,1)
        tint = psim.seeds(sidx,1);
        if tint~=tcur
            [tt,yt] = ode15s(ddt,[tcur,tint],yout(end,:)',opts);
            yend = yt(end,:)';
            tend = tt(end);
            yt([end],:) = interp1(tt,yt,[tt(end)-eps(tt(end))]);
            tt(end) = tt(end)-eps(tt(end));
            yout = [yout;yt];
            tout = [tout;tt];
        else
            yend = yout(:,end)';
            tend = tout(end);
        end
        y = yend;
        y(1) = y(1)+psim.seeds(sidx,2);        
        if tint~=tcur
            yout = [yout;y'];
            tout = [tout;tint];
        else
            yout(end) = y';
        end
        tcur = tint;
        sidx = sidx+1;
    elseif psim.purges(pidx,1)<=psim.seeds(sidx,1)
        tint = psim.purges(pidx,1);
        if tint~=tcur
            [tt,yt] = ode15s(ddt,[tcur,tint],yout(end,:)',opts);
            yend = yt(end,:)';
            tend = tt(end);
            yt([end],:) = interp1(tt,yt,[tt(end)-eps(tt(end))]);
            tt(end) = tt(end)-eps(tt(end));
            yout = [yout;yt];
            tout = [tout;tt];
        else
            yend = yout(:,end)';
            tend = tout(end);
        end
        y = yend;
        y(1) = y(1)-(psim.purges(pidx,2))*(y(1)/(y(1)+y(2)))*((1/ps.SVI)/(1-ps.watfrac));
        y(2) = y(2)-(psim.purges(pidx,2))*(y(2)/(y(1)+y(2)))*((1/ps.SVI)/(1-ps.watfrac));
        if tint~=tcur
            yout = [yout;y'];
            tout = [tout;tint];
        else
            yout(end) = y';
        end
        tcur = tint;
        pidx = pidx+1;
    else
        tint = psim.seeds(sidx,1);
        if tint~=tcur
            [tt,yt] = ode15s(ddt,[tcur,tint],yout(end,:)',opts);
            yend = yt(end,:)';
            tend = tt(end);
            yt([end],:) = interp1(tt,yt,[tt(end)-eps(tt(end))]);
            tt(end) = tt(end)-eps(tt(end));
            yout = [yout;yt];
            tout = [tout;tt];
        else
            yend = yout(:,end)';
            tend = tout(end);
        end
        y = yend;
        y(1) = y(1)+psim.seeds(sidx,2);        
        if tint~=tcur
            yout = [yout;y'];
            tout = [tout;tint];
        else
            yout(end) = y';
        end
        tcur = tint;
        sidx = sidx+1;
    end
end 
end

function dv = system_derivative(t,v,pb,ps,psim)
dv = zeros(9,1);
krep = psim.flow_fun(t)/ps.V_aer;
dv(1) = v(1)*v(3)*ps.k_gr_eng;
dv(2) = v(2)*v(3)*ps.k_gr_nat;
muec = ps.k_gr_eng*v(3);
munc = ps.k_gr_nat*v(3);
mur = -(dv(1)+dv(2))/(ps.V_aer-v(1)-v(2));
dv(3) = -v(3)*mur+...
    (-v(3))*(v(1)*ps.k_gr_eng+v(2)*ps.k_gr_nat)/(1-ps.watfrac)*1e6/ps.V_aer+...
    (psim.bod_in_fun(t)-v(3))*krep;
dv(4) = -v(4)*mur+...
    (psim.cr_in_fun(t)-v(4))*krep+...
    (-(ps.k_cr_imp_eng*v(1)+ps.k_cr_imp_nat*v(2))*v(4));
dv(5) = -v(5)*muec+...
    (ps.k_cr_imp_eng*v(4)*(ps.V_aer-v(1)-v(2)))+...
    -(v(5)*(pb.k_cat_cr_ox/3600)*v(9)*pb.m_crred/1000/(pb.K_M_cr_ox+v(5)));
dv(6) = -v(6)*munc+...
    (ps.k_cr_imp_nat*v(4)*(ps.V_aer-v(1)-v(2)));
dv(7) = -v(7)*muec+...
    pb.a_max_chr*1e6/(6.02e23)/(ps.v_cell/1e-15);
dv(8) = -v(8)*muec+...
    pb.a_max_trc2*(pb.K_M_lacI_dna/(pb.K_M_lacI_dna+v(7)))*1e6/(6.02e23)/(ps.v_cell/1e-15);
dv(9) = -v(9)*muec+...
    pb.a_max_chr*1e6/(6.02e23)/(ps.v_cell/1e-15);
end

function psim = gen_sim_params(varargin)
%   tspan is [t_start, t_end] in (hr)
%   flowfun (L/hr)
%   cr_in_fun (uM)
%   bod_in_fun (mg O2/L)
%   purges (purge table) (nx2 matrix) [times of purges (hr),
%       volumes of purges (L)] each row
fields = varargin(1:2:end);
values = varargin(1:2:end);
psim = struct('tspan',[0,24],'flow_fun',@(t)(25e6/24*3.78),...
    'bod_in_fun',@(t)(140),'purges',[],'cr_in_fun',@(t)(0.192),'seeds',[]);
if mod(nargin,2)~=0
    fprintf('Error: gen_sim_params_base must have even number of inputs')
end
for i = 1:(nargin/2)
    psim = setfield(psim,varargin{2*i-1},varargin{2*i});
end
end

function pb = gen_bioc_params(varargin)
pb = struct('k_cat_cr_ox',8812e-3*60,'K_M_cr_ox',41,'a_max_trc2',100,...
    'a_max_chr',100,'a_max_con',100,'n_p_mc',5,...
    'n_p_hc',30,'k_death',0.1,...
    'K_M_lacI_dna',0.18e-11,'K_M_chr_dna',0.18e-11,'m_crred',20400);
%copy numbers from https://openwetware.org/wiki/CH391L/S12/Origins_of_Replication
%lacI binding from http://pubs.acs.org/doi/pdf/10.1021/bi000924z
%ChrR mass from http://www.uniprot.org/uniprot/P0AGE6
%k_death arb
%assumed ChrR eq in spec to Lac
%all copy #'s max out at 100/hr
n = nargin;
if mod(n,2)~=0
    fprintf('Error: gen_bioc_params_base must have even number of inputs')
end
for i = 1:(n/2)
    pb = setfield(pb,varargin{2*i-1},varargin{2*i});
end
%{
pb.k_cat_cr_ox = %(umol/s)/(mg)(k_cat for cr oxidation)
pb.K_M_cr_ox = K_M_cr_ox;%(umol) (K_M for cr oxidation)
pb.a_max_trc2 = a_max_trc2;%(cp/hr)
pb.a_max_chr = a_max_chr;%(cp/hr)
pb.a_max_con = a_max_con;%(cp/hr)
pb.a_max_gntK = a_max_gntK;%(cp/hr)
pb.n_p_mc = n_p_mc;%(cp/cell)
pb.n_p_hc = n_p_hc;%(cp/cell)
pb.r_cell = r_cell;%cell radius in (um)
pb.l_cell = l_cell;%cell length in (um)
pb.k_death = k_death;%((fL)/(cp BamHI*hr))
pb.K_M_lacI_iptg = K_M_lacI_iptg;%(uM)
pb.h_lacI_iptg = h_lacI_iptg;%(dimensionless Hill exponent)
pb.K_M_lacI_dna = K_M_lacI_dna;%(uM)
pb.K_M_chr_cr = K_M_chr_cr;%(uM)
pb.h_chr_cr = h_chr_cr;%(dimensionless Hill exponent)
pb.K_M_chr_dna = K_M_chr_dna;(uM)
%}
end

function ps = gen_system_params(varargin)
ps = struct('k_gr_eng',1/12*0.5./24,'k_gr_nat',1/12/24,'V_aer',18.9e6,...
    'Y_xs',0.6e-6,'k_cr_imp_eng',0.10,'k_cr_imp_nat',0.01,...
    'watfrac',0.682,'SVI',108,'v_cell',1);
%assumes engineered bacteria grow at half-speed
%k_gr_nat from https://engineering.dartmouth.edu/~d30345d/courses/engs37/ActivatedSludge.pdf
%assumes Cr transport constant of 1 uM/hr, and 10x that in engineered cells
n = nargin;
if mod(n,2)~=0
    fprintf('Error: gen_system_params_base must have even number of inputs')
end
for i = 1:(n/2)
    ps = setfield(ps,varargin{2*i-1},varargin{2*i});
end
%{
%Y_xs from http://bionumbers.hms.harvard.edu/bionumber.aspx?id=105317&ver=4
ps.k_gr_eng = k_gr_eng;%growth in (hr^-1)/(mg O_2/L)
ps.k_gr_nat = k_gr_nat;%growth in (hr^-1)/(mg O_2/L)
ps.V_aer = V_aer;%volume of aerobic chamber (L)
ps.V_sed = V_sed;%volume of sedimentation chamber (L)
ps.E_sed_fun = E_sed_fun;%efficiency of sedimentation (C_i-C_o)/C_i (unitless)
ps.Y_xs = Y_xs;%mass-mass yield (mg dry biomass/mg BOD) (unitless)
ps.k_cr_imp_eng = k_cr_imp_eng;%efficiency of engineered CrVI uptake
    (umol/(L cell vol)/hr)
ps.k_cr_imp_nat = k_cr_imp_nat;%efficiency of natural  CrVI uptake
    (umol/(L cell vol)/hr)
ps.watfrac = watfrac;%proportion of sludge cells b.m. made of water http://pubmedcentralcanada.ca/pmcc/articles/PMC241608/pdf/aem00155-0081.pdf
ps.SVI = SVI;%v/m ratio of sludge dry weight occupancy of volume
%}
end
