function [Nd cref nref] = DynamicsMMSS_WM_NetworkDpMM(Nr,r0,rp,rint,CSD,A,B,KMM,ExtTh,DilTh,tauf,dtau)
%% Returns growth dynamics for communities of multiple species, based on the mechanistic model

% This code is free to use for any purposes, provided any publications resulting from the use of this code reference the original code/author.
% Author:  Babak Momeni (bmomeni@gmail.com); August 2016
% The code is not guaranteed to be free of errors. Please notify the author of any bugs, and contribute any modifications or bug fixes back to the original author.

%% Well-mixed model for growth of interacting species
% UIC: Uniform initial condition
% Ex: explicitly including the mediators
% MT: multi-target mediators
% rndseed = 1389;
% rand('twister',rndseed)
% DpMM: uptake based on Michaelis-Menten kinetics

% Nc = 2; % # of cell types
% Nm = 2; % # of mediators
% Nr = 300; % number of rounds of propagation
% r0 = [0.08; 0.1]; % population reproduction rates, per hour
% CSD = 0.9e8; % total initial cells
% Ksat = 1e7; % interaction strength saturation level of each population, fmole/ml
% K = 1e5; % Michaelis-Menten coefficient, fmole/ml 
% ExtTh = 0.1; % population extinction threshold
% DilTh = 1e8; % coculture dilution threshold
tau0 = 0;
% tauf = 250; % in hours
% dtau = 0.01; % in hours, cell growth update and uptake timescale
% at = 0.1; % avg. consumption values (fmole per cell); alpha_ij: population i, resource j
% bt = 0.1; % avg. production rates (fmole per cell per hour); beta_ij: population i, resource j

% rint = [0 0.1; 0.05 0]; % matrix of interaction coefficients
[Nc Nm] = size(rint);
% Si = Ksat*(0.5+rand(Nm,1)); % interaction saturation concentration, fmole/ml
% KMM = K*(0.5+rand(Nm,Nc)); % Michaelis-Menten coefficients, fmole/ml

%% Parameters
% R = [0 1; 1 0];
% P = [1 0; 0 1];

% interaction matrix
% alpha = at*(0.5+rand(Nc,Nm)); % consumption rates
% beta = bt*(0.5+rand(Nc,Nm)); % mediator release rates
% A = (R.*alpha)';
% B = (P.*beta)';

%% Initial state 
% rp = 1/Nc*ones(1,Nc); % cell distrbution
C = zeros(Nm,1); % concentrations of interaction mediators, fmole/ml

%% Cell-growth time-course
taurng = tau0:dtau:tauf;
Nt = length(taurng);
Nd = zeros(1,Nr);

tc = zeros(1,Nr*Nt);
Xc = zeros(Nc,Nr*Nt);
Cc = zeros(Nm,Nr*Nt);

cct = 0;
X = CSD*rp'; % initial number of each cell type
for rnd = 1:Nr,
    C = CSD/sum(X)*C;
    X = CSD*rp'; % initial number of each cell type
    
    tau0 = 0; % in hours
    tau = tau0;
    
    Xlm = zeros(Nc,Nt);
    Clm = zeros(Nm,Nt);
    rzm = zeros(Nc,Nt);
    ct = 0;
    while (tau<=tauf-dtau)&&(sum(Xlm(:,max(ct,1)))<DilTh),
        
        ct = ct+1;
        tau = taurng(ct);

        Ce = (C.*(C>0))*ones(1,Nc);
        rz = r0 + (rint./(Ce+KMM)')*(C.*(C>0));
        X = X + dtau*(rz.*X);
        X(X<ExtTh) = 0;
        AMM = A'.*Ce./(Ce+KMM);
        C = C + dtau*(B'*X - AMM*X);
        % C(C<0) = 0;
        Xlm(:,ct) = X;
        Clm(:,ct) = C(:,1);
        rzm(:,ct) = rz;
        
    end
    Nd(rnd) = ct;
    rp = 1/sum(Xlm(:,ct))*Xlm(:,ct)';
    
    if cct ==0,
        tc(cct+1:cct+ct) = taurng(1:ct);
    else
        tc(cct+1:cct+ct) = tc(cct) + taurng(1:ct);
    end
    Xc(:,cct+1:cct+ct) = Xlm(:,1:ct);
    Cc(:,cct+1:cct+ct) = Clm(:,1:ct);
    cct = cct+ct;

end
cref = Cc(1:Nm,1:cct);
nref = Xc(1:Nc,1:cct);

% figure
% semilogy(tc,Xc')
% % hold on
% % semilogy(Xlm(Ne,:)','o')
% % xlim([0 1000])
% xlabel('Time (hrs)')
% ylabel('Species abundance (cells)')
% set(gca,'YTick',logspace(-4,10,8))
% 
% figure
% semilogy(tc,1e-6*Cc')
% % hold on
% % semilogy(Clm(Ne,:)','o')
% % xlim([0 1000])
% ylim([1e-6 1e6])
% xlabel('Time (hrs)')
% ylabel('Chemical concentrations (\muM)')
% set(gca,'YTick',logspace(-6,10,9))

return;