Competing Yeast Species

Kurt Bryan and SIMIODE

A notebook 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 "nd" for "no data").

sacc1 = {0.37, 8.87, 10.66, 12.5, 13.27, 12.87, 12.70, nd, nd, nd, nd}

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 = {nd, 1, nd, 1.7, nd, 2.73, nd, 4.87, 5.67, 5.8, 5.83}

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

schiz2 = {nd, 1.27, nd, 2.33, nd, nd, 4.56}

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

saccmix1 = {0.375, 3.99, 4.69, 6.15, nd, 7.27, 8.3, nd, nd, nd, nd}

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, nd, 1.71, 1.84, nd, nd, nd, nd}

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 for runs 1 and 2, for the Saccharomyces population alone (column 2), the Schizosaccharomyces population alone (column 4), the Saccharomyces population in competition (column 3), and the Schizosaccharomyces population in competition (column 5).

The data will be in (time, population) pairs.

N1 = Length[times1]; N2 = Length[times2];

saccind = Join[Table[{times1[[j]],sacc1[[j]]},{j,1,N1}],Table[{times2[[j]],sacc2[[j]]},{j,1,N2}]];

schizind = Join[Table[{times1[[j]],schiz1[[j]]},{j,1,N1}],Table[{times2[[j]],schiz2[[j]]},{j,1,N2}]];

saccmix = Join[Table[{times1[[j]],saccmix1[[j]]},{j,1,N1}],Table[{times2[[j]],saccmix2[[j]]},{j,1,N2}]];

schizind = Join[Table[{times1[[j]],schizmix1[[j]]},{j,1,N1}],Table[{times2[[j]],schizmix2[[j]]},{j,1,N2}]];

We can plot the data. For example, for the saccharomyces population growing in isolation, with all "nd" points replaced by -1 and the view restricted to nonnegative population values so those points don't show up.

datplot = saccind/.{nd->-1};
ListPlot[datplot,PlotRange->{{0,60},{0,15}}]

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

logde = u'[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[[1,2]]

Solve the logistic ODE with this data and assign it to the function "u1sol":

sol = DSolve[{logde,u[t0]==u0},u,t];
u1sol=u/. sol[[1]]

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

SS = 0;
N0 = N1 + N2;
For[i = 1, i ≤N0, i++,
If[ NumericQ[saccind[[i,2]]], (* Check to see that saccind[[i,2]] is defined *)
SS = SS +( u1sol[saccind[[i,1]]]-saccind[[i,2]])^2
];
]

Minimize the sum of squares to obtain a best fit for r1 and K1, using Mathematica's  FindMinimum command. Initial guess is K1 = 12, r1 = 0.5.

parmin=FindMinimum[SS,{{K1,12},{r1,0.5}}]

Define the best-fit solution:

bestu1[t_] = u1sol[t]/.parmin[[2]]

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

datplot = saccind/.{nd->-1};
plt1 = ListPlot[datplot,PlotRange->{{0,60},{0,15}}];
plt2 = Plot[bestu1[t],{t,0,60},PlotStyle->Red];
Show[plt1,plt2]

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

de1 = u1'[t] == r1*u1[t]*(1-u1[t]/K1-a*u2[t]/K1)
de2 = u2'[t] == r1*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.