%% Competing Yeast Species
%% Kurt Bryan and SIMIODE
% A script to facilitate analyzing the competing yeast species project in Section 7.5.3 of the text. The goal is to use the data of Table 7.2 to estimate the parameters r1, K1, r2, K2, a, and b in equations (7.3)-(7.4).
% The Data:
% Time (hours) for yeast populations, run 1 (first 11 lines of Table 7.2).

times1 = [6, 16, 24, 29, 40, 48, 53, 72, 93, 117, 141];

% And for run 2 (remaining 7 lines in Table 7.2)

times2 = [7.5, 15, 24, 31.5, 33, 44, 51.5];

% Saccharomyces population, run 1 (column 2 in Table 7.2, with "-1" indicating  "no data").

sacc1 = [0.37, 8.87, 10.66, 12.5, 13.27, 12.87, 12.70, -1, -1, -1, -1];

% Saccharomyces population, run 2, (column 2, bottom 7 lines in Table 7.2)

sacc2 = [1.63, 6.2, 10.97, 12.60, 12.9, 12.77, 12.9];

% Schizosaccharomyces population run 1, (column 4 in Table 7.2, with "nd" for "no data").

schiz1 = [-1, 1, -1, 1.7, -1, 2.73, -1, 4.87, 5.67, 5.8, 5.83];

% Schizosaccharomyces population run 2, (column 4, bottom 7 lines in Table 7.2)

schiz2 = [-1, 1.27, -1, 2.33, -1, -1, 4.56];

% Saccharomyces population, run 1, mixed Sach/Schiz experiment (column 3, first 11 lines)

saccmix1 = [0.375, 3.99, 4.69, 6.15, -1, 7.27, 8.3, -1, -1, -1, -1];

% Saccharomyces population, run 2, mixed Sach/Schiz experiment (column 3, bottom 7 lines)

saccmix2 = [0.923, 3.082, 5.78, 9.91, 9.47, 10.57, 9.883];

% Schizosaccharomyces population run 1, mixed Sach/Schiz experiment (column 5, first 11 lines)

schizmix1 = [0.291, 0.98, 1.47, 1.46, -1, 1.71, 1.84, -1, -1, -1, -1];

% Schizosaccharomyces population run 2, mixed Sach/Schiz experiment (column 5, bottom 7 lines)

schizmix2 = [0.371, 0.63, 1.22, 1.112, 1.225, 1.102, 0.961];

% Let's form amalgamated data sets, by putting the runs together in (time, population) pairs

N1 = length(times1); N2 = length(times2); N0 = N1+N2;
times = [times1 times2];
saccind = zeros(2,N0);
saccind(1,:)=times; saccind(2,1:N1)=sacc1;  saccind(2,N1+1:N0)=sacc2;
schizind = zeros(2,N0);
schizind(1,:)=times; schizind(2,1:N1)=schiz1;  schizind(2,N1+1:N0)=schiz2;
saccmix = zeros(2,N0);
saccmix(1,:)=times; saccmix(2,1:N1)=saccmix1;  saccmix(2,N1+1:N0)=saccmix2;
schizmix = zeros(2,N0);
schizmix(1,:)=times; schizmix(2,1:N1)=schizmix1;  schizmix(2,N1+1:N0)=schizmix2;

% We can plot the saccharomyces data with

scatter(saccind(1,:),saccind(2,:),'filled')
axis([0 60 0 15])

% The restriction of the axes prevents the "no data" -1 points from showing up.
% A similar plot can be constructed for the schizosaccharomyces population in isolation, or either species in competition.
% Estimate Growth Parameters for Saccharomyces: We can estimate the parameters r1 and K1 for the logistic growth model for the saccharomyces population using the corresponding data.
% The logistic ODE is

syms u(t);
syms r1 K1;
logde = diff(u(t),t) == r1*u(t)*(1-u(t)/K1)

% The initial data could be taken as any (time, population) data point, but we'll use the first:

t0 = saccind(1,1);
u0 = saccind(2,1);

% Solve the logistic ODE with this data

u1sol(t) = dsolve(logde,u(t0)==u0)

% A plot or even a glance at the data suggests that K1 is around 12 to 14. To estimate K1 and r1 more carefully, form a sum of squares that ignores missing data.

syms SS;
SS = 0;
for j=1:N0
    if saccind(2,j)>0
        SS = SS + (u1sol(saccind(1,j))-saccind(2,j))^2;
    end
end

% Minimize the sum of squares to obtain a best fit for r1 and K1, using Matlab's "fminsearch" command. Initial guess is K1 = 12, r1 = 0.5. We will use Matlab's fminsearch algorithm to seek the minimizer, after converting "SS" to a Matlab function "SSm":

SSm = matlabFunction(SS,'Vars',{[r1 K1]});
[optpars,fval] = fminsearch(SSm,[0.5 12.0])

% Define the best-fit solution:

u1best(t) = subs(u1sol(t),[r1 K1],optpars);

% Compare the plot and the data on an appropriate range (up to t = 60 here).

scatter(saccind(1,:),saccind(2,:),'filled')
axis([0 60 0 15])
hold on
fplot(u1best(t),[0 60])
hold off

% A similar procedure can be used to estimate r2 and K2 for the schizosaccharomyces data, using (for example) the data point t0 = 16, u0 = 1 as an "initial" data point.
% The Competition Parameters a and b: The ODEs that govern the competing species system are

syms u1(t) u2(t);
syms r2 K2 a b;
de1 = diff(u1(t),t) == r1*u1(t)*(1-u1(t)/K1-a*u2(t)/K1)
de2 = diff(u2(t),t) == r2*u2(t)*(1-u2(t)/K2-b*u1(t)/K2)

% A closed-form solution cannot be obtained, but a guess-and-plot approach based on solving these ODEs numerically can work, using the values for r1, K1, r2, K2 found above.