Two Phytoplankton Two Resource Model

% two resource, two phytoplankton, chemostat
clear all; close all;

% initialize figure plotting...
figure
set(gcf,'Renderer','zbuffer') 


% timestep (days)
deltat = 0.1;
% number of days to timestep
max_days = 150.0 ;
nstepmax = max_days/deltat;
% nout = how often to plot data point (integer days)
nout = 0.5/deltat ;

% set phyto characteristics
% growth (day-1)
  mu1 = 0.50;
  mu2 = 0.50;
% diluation rate (day-1)
  dil =  0.2;
% half-saturation, K (micromol l-1)
  k1a =  0.1;
  k1b =  0.2;
  k2a =  0.2;
  k2b =  0.1;

% elemental ratios -  resource b:resource a required by phyto
 R1ba = 2.0 ;
 R2ba = 1.0 ;

% incoming resource concentration (micromol l-1)
  RINba = 1.7 ;
  nutina = 0.5 ;
  nutinb = RINba*nutina ;

% initialization 
% resources (micromol l-1)
  nuta = nutina;
  nutb = nutinb;
% phytoplankton (micromol l-1)
  phy1 = 0.002;
  phy2 = 0.002;
% initial time = 0.0
tdays = 0.0;

% begin timestepping.....................
for nstep = 1:nstepmax,
  tdays = tdays + deltat;

% determine tendencies

% net nutrient supply rate (micromol l-1 day-1) 
  sna = dil*(nutina - nuta) ;
  snb = dil*(nutinb - nutb) ;

% phytoplankton 1
  limit_1a = nuta/(nuta+k1a);
  limit_1b = nutb/(nutb+k1b);  
  limit1 = min(limit_1a, limit_1b) ;
  grow1 = mu1*limit1*phy1;
  loss1 = dil*phy1;
  dphy1dt = grow1 - loss1 ;

% phytoplankton 2 
  limit_2a = nuta/(nuta+k2a);
  limit_2b = nutb/(nutb+k2b);  
  limit2 = min(limit_2a, limit_2b) ;
  grow2 = mu2*limit2*phy2;
  loss2 = dil*phy2;
  dphy2dt = grow2 - loss2 ;

% resources
   dnutadt = - (grow1+grow2) + sna;
   dnutbdt = - (grow1*R1ba + grow2*R2ba) + snb;

% forward step
  phy1 = phy1 + dphy1dt*deltat;
  phy2 = phy2 + dphy2dt*deltat;
  nuta = nuta + dnutadt*deltat;
  nutb = nutb + dnutbdt*deltat;


% update variables for plotting
% only plot at interval defined by nstep
if mod(nstep, nout) == 0
% phytoplankton
  phystore1(nstep) = phy1;
  phystore2(nstep) = phy2;
% resources 
  nutstorea(nstep) = nuta;
  nutstoreb(nstep) = nutb;
% time
  tstore(nstep) = tdays;
end

end
% end timestepping ......................


% plot results 
    figure(gcf);
    plot(tstore,phystore1,'.b');
    set(gca,'FontSize',14) ;
    axis([0 tdays 0.0 1.5]);
    xlabel('time (days) ','FontSize',16);
    ylabel('concentration  (umol l^-^1) ','FontSize',16);
    grid;
    hold on;
    plot(tstore,phystore2,'.k');
    plot(tstore,nutstorea,'.r');
    plot(tstore,nutstoreb,'.g');
    hleg1 = legend('phyto1','phyto2','nutrient a','nutrient b');