function nLV = DynamicsWM_NetworkBFSatLV(Nr,rpLV,rint_SatLV,CSD,Si_SatLV,ExtTh,Nd,dtau,taurng,r0_est,rs_est,Si_est)
%% Returns growth dynamics for communities of multiple species, based on the canonical pairwise 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.

Nc = length(rpLV);
Nt = sum(Nd);
tcLV = zeros(1,Nr*Nt);
XLVc = zeros(Nc,Nr*Nt);
cct = 0;
for rnd = 1:Nr,
    XLV = CSD*rpLV'; % initial number of each cell type
    
    XLVlm = zeros(Nc,Nt);
    ctLV = 0;
    while (ctLV<Nd(rnd)),
        ctLV = ctLV+1;

        rLVz = diag(r0_est + rs_est.*XLV./(XLV+Si_est)) + (rint_SatLV./(ones(Nc,1)*XLV'+Si_SatLV)).*(XLV*ones(1,Nc));
        XLV = XLV + dtau*(rLVz*XLV);
        XLV(XLV<ExtTh) = 0;
        XLVlm(:,ctLV) = XLV;
%         if mod(ctLV,1000)==1,
%             disp(rLVz)
%         end
    end
    rpLV = 1/sum(XLVlm(:,ctLV))*XLVlm(:,ctLV)';
%     disp(rpLV)
    if cct ==0,
        tcLV(cct+1:cct+ctLV) = taurng(1:ctLV);
    else
        tcLV(cct+1:cct+ctLV) = tcLV(cct) + taurng(1:ctLV);
    end
    XLVc(:,cct+1:cct+ctLV) = XLVlm(:,1:ctLV);
    cct = cct+ctLV;

end
nLV = XLVc(1:Nc,1:cct);

% figure
% semilogy(taurng(1:cct),n1LV(1:cct),'b:')
% hold on
% semilogy(taurng(1:cct),n2LV(1:cct),'r:')
% xlabel('Time (hrs)')
% ylabel('Population size')
% ylim([1e2 1e10])

