% 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 rint A B KMM
rndseed = 4000;
rand('twister',rndseed)

r0 = [0.12; 0.1]; % population reproduction rates, per hour
at = 0.05; % avg. consumption values (fmole per cell); alpha_ij: population i, resource j
bt = 0.5; % avg. production rates (fmole per cell per hour); beta_ij: population i, resource j
rint = [0; 0.03]; % Nc*Nm matrix of interaction coefficients

CSD = 1e5; % total initial cells
K = 1e1; % 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
Nr = 10; % number of rounds of propagation

[Nc Nm] = size(rint);
KMM = K*[1 1]; % Michaelis-Menten coefficients for consumption and influence, fmole/ml

rp0 = [0.5 0.5]; %1/Nc*ones(1,Nc);

%% Connectivity, Nm*Nc
R = [0 1]; % release
P = [1 0]; % consumption

% interaction matrix
alpha = at*[0 1]; % consumption rates
beta = bt*[1 0]; % mediator release rates
A = (R.*alpha)';
B = (P.*beta)';

%% Initial state
rp = rp0; % cell distrbution

[tout,yout] = ode23(@Dynamics_WM_NetworkDpMM_ODE23,[0 58],[rp(1)*CSD rp(2)*CSD 0]);
nref = yout(:,1:2)';
cref = yout(:,3)';

nrefs = zeros(Nr,600000);
crefs = zeros(Nr,600000);
cct = 0;
cf = 1;
Nd = zeros(1,Nr);
for nr = 1:Nr-1,
    disp(nr)
    tDil = interp1(sum(nref)',tout,DilTh);
    Nd(nr) = sum(tout<tDil);
    nrefc(:,cct+1:cct+Nd(nr)) = cf*nref(:,1:Nd(nr));
    crefc(:,cct+1:cct+Nd(nr)) = cref(:,1:Nd(nr));
    nrefs(nr,1:Nd(nr)) = cf*nref(:,1:Nd(nr));
    crefs(nr,1:Nd(nr)) = cref(:,1:Nd(nr));
    cf = cf*sum(nref(:,Nd(nr)))/CSD;
    cct = cct + Nd(nr);
    rp = 1/sum(nref(:,Nd(nr)))*nref(:,Nd(nr))';
    C0 = cref(Nd(nr))*CSD/sum(nref(:,Nd(nr)));
    [tout,yout] = ode23(@Dynamics_WM_NetworkDpMM_ODE23,[0 58],[rp(1)*CSD rp(2)*CSD C0]);
    nref = yout(:,1:2)';
    cref = yout(:,3)';
end

tDil = interp1(sum(nref)',tout,DilTh);
Nd(Nr) = sum(tout<tDil);
nrefc(:,cct+1:cct+Nd(Nr)) = cf*nref(:,1:Nd(Nr));
crefc(:,cct+1:cct+Nd(Nr)) = cref(:,1:Nd(Nr));
nrefs(Nr,1:Nd(Nr)) = cf*nref(:,1:Nd(Nr));
crefs(Nr,1:Nd(Nr)) = cref(:,1:Nd(Nr));

figure
semilogy(tout,nref)
xlabel('Time (hrs)')
ylabel('Population densities')
% xlim([0 Ng])
% ylim([1e-2 1.2])
% set(gca,'YTick',[1e-2 1e-1 1])

figure
semilogy(1/log(2)*(log(sum(nref,1))-log(CSD/Nc)),1./(ones(Nc,1)*sum(nref,1)).*nref)
xlabel('Community doublings')
ylabel('Population fraction')
% xlim([0 Ng])
ylim([1e-2 1.2])
set(gca,'YTick',[1e-2 1e-1 1])

taurngc = 0:dtau:(cct-1)*dtau;
figure
semilogy(1/log(2)*(log(sum(nrefc,1))-log(CSD/Nc)),1./(ones(Nc,1)*sum(nrefc,1)).*nrefc)
xlabel('Community doublings')
ylabel('Population fraction')
% xlim([0 Ngo])
ylim([1e-2 1.2])
% set(gca,'YTick',[1e-3 1e-2 1e-1 1])

figure
semilogy(1/log(2)*(log(sum(nrefc,1))-log(CSD/Nc)),crefc)
hold on
semilogy(smooth(1/log(2)*(log(sum(nrefc,1))-log(CSD/Nc)),100),smooth(crefc,100),'r')
xlabel('Community doublings')
ylabel('Mediator concentration')
xlim([0 100])
ylim([1e-2 1e2])
% set(gca,'YTick',[1e-3 1e-2 1e-1 1])

save(strcat('CalculateDynamicsODE23_C1Sp2_CmnsDp_0P_a',num2str(at),'_b',num2str(bt),'_rp_',num2str(rp0(1)),'_rint',num2str(rint(2,1)),'_CSD',num2str(CSD),'_DilTh',num2str(DilTh),'_Nr',num2str(Nr),'_rndseed',num2str(rndseed),'.mat'))
