function [rint_SatLV Si_SatLV] = DeriveBFSatLVMMSS_WM_DpMM_r21(Nr,r0,rint,CSD,A,B,KMM,tauf,r0_est,rs_est,Si_est)
%% Estimates canonical pairwise model interaction parameters (Step 2)

% 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.

%% Deriving equivalent fitness model for chemically interacting species
% r0: population reproduction rates, per hour
% rint: matrix of interaction coefficients
% A: consumption matrix
% B: production matrix
% Si: population levels for influence saturation
% CSD: total initial cells
% dtau: simulation time-step
% tauf: total simulation time

global n1 n2 Nt taufe dtau ExtTh DilTh r0i rsi Sii


% Nc: # of cell types
% Nm: # of mediators
Nc = size(A,1);

%% Cell-growth time-course
tau0 = 0;
taurng = tau0:dtau:tauf;

rint_SatLV = zeros(Nc,Nc);
Si_SatLV = zeros(Nc,Nc);
for ii = 1:Nc,
    for jj = ii+1:Nc,
        % focal pair

        %% Initial state
        rpp = zeros(1,Nc);
        rpp([ii jj]) = [0.5 0.5]; % initial composition

        [Nd cref nref] = DynamicsMMSS_WM_NetworkDpMM(Nr,r0,rpp,rint,CSD,A,B,KMM,ExtTh,DilTh,tauf,dtau);
        n1 = nref(ii,sum(Nd(1:Nr-1))+1:sum(Nd(1:Nr)));
        n2 = nref(jj,sum(Nd(1:Nr-1))+1:sum(Nd(1:Nr)));
        taufe = taurng(Nd(Nr));

%         figure
%         semilogy(taurng(1:Nd(Nr)),n1,'b')
%         hold on
%         semilogy(taurng(1:Nd(Nr)),n2,'r')
%         xlabel('Time (hrs)')
%         ylabel('Populations')
%         pause
        
        x0 = [4 4 6 6];
        r0i = r0_est([ii jj]);
        rsi = rs_est([ii jj]);
        Sii = Si_est([ii jj]);
        Nt = Nd(Nr);
        x = lsqnonlin(@FitCost_BFSatLV_r21,x0);

        rint_SatLV(ii,jj) = 0.01*x(1);
        rint_SatLV(jj,ii) = 0.01*x(2);
        Si_SatLV(ii,jj) = 10.^x(3);
        Si_SatLV(jj,ii) = 10.^x(4);
%         disp([ii jj x])
    end
end

return;
