% Fit dynamics of ExMT model assuming SatLV model to compare with DpMM

% SatLV: saturating Lotka-Volterra dr_i = phi_ij * n_j/(n_j + K_j)
% DpMM: mechanistic with consumable Michaelis-Menten interaction mediators
% FitDynamicsMM: Same Michaelis-Menten coefficients used both for consumption and for influence

clear

global r0 dtau ExtTh DilTh CSD rp0 n1 n2 Nt taufe rii Kii
rndseed = 2611;
rand('twister',rndseed)

r0 = [-0.02; -0.01]; % population reproduction rates, per hour
CSD = 10; % total initial cells
K = 1e5; % Michaelis-Menten coefficient, fmole/ml 
ExtTh = 1e-16;
DilTh = 1e16; % coculture dilution threshold
tau0 = 0;
tauf = 5000; % in hours
dtau = 0.01; % in hours, cell growth update and uptake timescale
at = 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
Nr = 1; % number of rounds of propagation
KC = K*[1 0;
    2 1]; % Nc*Nm 
 
rc = [15; 25]; % Nc*1 matrix of resource-mediated fitness benefit
Nc = length(r0);
Nm = length(KC);
 
rp0 = 1/Nc*ones(1,Nc);
% rp0 = [1 0];
% rp0 = [0 1];
 
%% Connectivity, Nm*Nc
P = [1 1; 0 1]; % consumption
R = [0 0; 1 0]; % release
%% Rates, Nm*Nc
alpha = at*P; % consumption rates
beta = bt*[0 0; 10 0]; % mediator release rates
beta0 = bt*[10 0]; % resource supply rates

%% interaction matrix, Nc*Nm
A = (P.*alpha)'; % consumption matrix
B = (R.*beta)'; % release matrix

%% Initial state 
rp = rp0; % cell distrbution

[Nd cref nref] = Dynamics_WM_NetworkDp_ARN(Nr,r0,rp,rc,CSD,A,B,beta0,KC,DilTh,tauf,dtau);

% figure
% semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nref(1,1:sum(Nd(1:Nr)))./sum(nref(:,1:sum(Nd(1:Nr)))),'b-')
% hold on
% semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nref(2,1:sum(Nd(1:Nr)))./sum(nref(:,1:sum(Nd(1:Nr)))),'-','color',[0 0.5 0])
% xlabel('Time (hrs)')
% ylabel('Population fraction')
% ylim([0.01 1.1])
% 
figure
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,cref(1,1:sum(Nd(1:Nr))),'k-')
hold on
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,cref(2,1:sum(Nd(1:Nr))),'-','color',[0.6 0.6 0.6])
ylim([1e-4 1e4])
set(gca,'YTick',[1e-4 1e-2 1 1e2 1e4])
xlabel('Time (hrs)')
ylabel('Mediator concentrations (\muM)')

rii = [0; 0]; % [r0(1,1)+rc(1,1); 0];
Kii = [-beta0(1,1)/alpha(1,1)/r0(1,1); 0];
rint_LogLV = [0 0; 0 0]; 

% rpLV = [0.5 0.5]; % cell distrbution
rpLV = rp0; % cell distrbution
taurng = 0:dtau:tauf;

xest0 = [log10(5e2) 0.05 0 1]; % initial guess, [log10(Kii(2,1)) 100*rint_LogLV(2,1)]
Nest = 150/dtau; % range of time points to be used for fitting
Nt = Nest;
n1 = nref(1,1:Nest);
n2 = nref(2,1:Nest);
taufe = taurng(Nest);
xest = lsqnonlin(@FitCost_BFLogLV_ARN_DI,xest0);
Kii(2,1) = 10^xest(1);
rint_LogLV(2,1) = 0.01*xest(2);
dii(1,1) = 0.01*xest(3);
dii(2,1) = 0.01*xest(4);

Kii
rint_LogLV
dii

nLV = DynamicsWM_NetworkBFLogLV_DI(Nr,rpLV,rint_LogLV,CSD,ExtTh,Nd,dtau,taurng,rii,Kii,dii);
nLV1 = nLV(1,:);
nLV2 = nLV(2,:);

figure
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nref(1,1:sum(Nd(1:Nr))),'b-')
hold on
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nref(2,1:sum(Nd(1:Nr))),'-','color',[0 0.5 0])
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nLV1,'b:')
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nLV2,':','color',[0 0.5 0])
xlabel('Time (hrs)')
ylabel('Population size')

figure
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nref(1,1:sum(Nd(1:Nr)))./sum(nref(:,1:sum(Nd(1:Nr))),1),'b-')
hold on
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nref(2,1:sum(Nd(1:Nr)))./sum(nref(:,1:sum(Nd(1:Nr))),1),'-','color',[0 0.5 0])
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nLV1./sum(nLV,1),'b:')
semilogy(0:dtau:(sum(Nd(1:Nr))-1)*dtau,nLV2./sum(nLV,1),':','color',[0 0.5 0])
ylim([1e-4 1.2])
xlabel('Time (hrs)')
ylabel('Population fraction')

% figure
% loglog(nref(1,1:sum(Nd(1:Nr))),nref(2,1:sum(Nd(1:Nr))),'-')
% xlabel('Population 1')
% ylabel('Population 2')

save('C2Sp2_ARCLi_NoSatDp_FitBFLogLV_DI_ConstSupply_Resurge3.mat')
