function LiverFit = ModelSimulation_KineticsThresholds
dose=2*10^6
function C=kinetics(theta,t)
f1=20
f2=10
f3=0.002
k1_param=0.02
k2=0.02
r1=0.196
r2=0.095
r3=0.110
c1=3.17
c2=2.59
c3=3.84
M1=10^9.55
M2=10^7.9
M3=10^4.4
gamma1=0.00074
gamma2=0.00074
gamma3=0.0014

c0=[(20/(20+10+0.002))*dose;(10/(20+10+0.002))*dose;(0.002/(20+10+0.002))*dose];
[T,Cv]=ode89(@DifEq,t,c0);
    function dC=DifEq(t,c)
    dcdt=zeros(3,1);
    dcdt(1)= (f1/(f1+f2+f3))*k1_param.*c(1)+(f1/(f1+f2+f3))*k2.*c(2)+(r1).*(c(1).*(1-(c(1)/(M1)))-k1_param.*c(1)-(c1).*c(1)/(1+gamma1*c(1)));
    dcdt(2)= (f2/(f1+f2+f3))*k1_param.*c(1)+(f2/(f1+f2+f3))*k2.*c(2)+(r2).*(c(2)*(1-(c(2)/(M2)))-k2.*c(2)-(c2).*c(2)/(1+gamma2*c(2)));
    dcdt(3)= (f3/(f1+f2+f3))*k1_param.*c(1)+(f3/(f1+f2+f3))*k2.*c(2)+(r3).*(c(3).*(1-(c(3)/(M3)))-(c3).*c(3)/(1+gamma3*c(3)));
    dC=dcdt;
    end
C=Cv;
end

%The following is the data obtained by Join-Lambert et al. used in fitting
%the model
t = [0,8,24,48];
Liver = [0;101329.8755;9983499.7496;3397799162.7540];
Spleen = [0;365778.6821;7088832.4672;169779424.6359];
Brain = [0;0;10^3.04;10^4.62];
c = [Liver, Spleen, Brain];

theta0=[1;1;1;1];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta);
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Fit = kinetics(theta, tv);
fig = figure;

LiverFit = Fit(:,1);
SpleenFit = Fit(:,2);
BrainFit = Fit(:,3);

subplot(3,1,1);
plot(tv, LiverFit, 'r');
hold on 
plot(t, Liver, 'bo')
title('Figure 3a')
lgd = legend('Liver (t) Prediction', 'Join-Lambert data', 'location', 'northwest')
lgd.FontSize = 6;

subplot(3,1,2);
plot(tv, SpleenFit,'r');
hold on 
plot(t, Spleen,'bo')
title('Figure 3b')
lgd = legend('Spleen(t) Prediction', 'Join-Lambert data', 'location', 'northwest')
lgd.FontSize = 6;

subplot(3,1,3);
plot(tv, BrainFit,'r');
hold on 
plot(t, Brain,'bo')
title('Figure 3c')
ylim([0 10^4.7])
lgd = legend('Brain(t) Prediction', 'Join-lambert data', 'location', 'northwest')
lgd.FontSize = 6;

han=axes(fig,'visible','off'); 
han.XLabel.Visible='on';
han.YLabel.Visible='on';
ylabel(han,'Bacteria in liver, spleen, and brain respectively (CFU)');
xlabel(han,'Time post infection (Hours)');
end