Two Phytoplankton One Resource Model

% one resource, two phytoplankton    

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.01/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)
  k1 =  0.1;
  k2 =  0.2;

 
% non-linear loss, "grazing" ( (micromol l-1 day)-1  )
  g1 = 0.4
  g2 = 0.2
  
% incoming resource concentration (micromol l-1)
  nutin1 = 1.0 ;

% initialization 
% resources (micromol l-1)
  nut1 = nutin1;
% 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
  n0 = 1.0;
  nutin1 = n0;
% sinusoidally varying resource supply
% changing concentration of inflowing medium
% UNCOMMENT LINE BELOW FOR VARYING NUTRIENT SOURCE
% nutin1 = n0*( 1.0 + 0.5*(sin(2.0*pi*tdays/5.0)) ); 
  
% net nutrient supply rate (micromol l-1 day-1) 
  sn1 = dil*(nutin1 - nut1) ;

% phytoplankton 
  grow1 = mu1*nut1/(nut1+k1)*phy1;
  loss1 = dil*phy1;
  graz1 = 0.0;
% UNCOMMENT LINE BELOW FOR QUADRATIC LOSS TERM
% graz1 = g1*phy1*phy1;
  dphy1dt = grow1 - loss1 - graz1 ;

  grow2 = mu2*nut1/(nut1+k2)*phy2;
  loss2 = dil*phy2;
  graz2 = 0.0;
% UNCOMMENT LINE BELOW FOR QUADRATIC LOSS TERM
%  graz2 = g2*phy2*phy2 ;
  dphy2dt = grow2 - loss2 - graz2;

% resources
   dnut1dt = - (grow1+grow2) + sn1;

% forward step
  phy1 = phy1 + dphy1dt*deltat;
  phy2 = phy2 + dphy2dt*deltat;
  nut1 = nut1 + dnut1dt*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 
  nutstore1(nstep) = nut1;
% time
  tstore(nstep) = tdays;
end

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


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