from collections import OrderedDict
import copy
import numpy as np
import matplotlib.pyplot as plt
import os
import scipy.optimize as spopt
import platereaderdataloader as prdl

%matplotlib inline

Table of Contents

• Description of the Evolution Experiment
• Evolution of the growth rate
• Evolution of saturating OD

Description of the Evolution Experiment

The evolution experiment (short compared to Lenski) was similar to Lenski's experiment. I started with a single strain of E. coli and split it into 45 populations that each grew in 500 microliters of supplemented M9 medium at 30 degrees Celsius in a platereader each day. There were also 3 blank wells. Over the course of the day, the E. coli consume their food and the platereader measures the optical density of the E. coli every ~10 minutes. At the end of each day, all 45 populations are diluted 500x and put into a fresh well of 500 microliters then back in the platereader.

Every three days, a sample of all populations is frozen in case of mishap. Twice the experiment was stopped, all stocks were frozen, and then the experiment was resumed later. Let's look at the growth curve data.

What makes this experiment different from Lenski's is that the strain of E. coli used has an adjustable mutation rate tuned by otherwise innocuous chemical inducers. 5 combinations of inducer concentrations were chosen giving 5 mutation rates with 9 replicates of each. Thus, for an otherwise fixed strain of E. coli we can see how the mutation rate affects evolution.

First I need to import the data.

os.chdir("Experiment 1")

days = OrderedDict()

filename_end = ' July 2018 NS001 delta CAT 30 C.xlsx'
platelayout1 = ['High '+str(i) for i in range(1,10)] + ['HiMid '+str(i) for i in range(1,10)] + \
               ['Mid '+str(i) for i in range(1,10)] + ['LoMid '+str(i) for i in range(1,10)] + \
               ['Low '+str(i) for i in range(1,10)] + ['blank1', 'blank2', 'blank3']

for i in range(9,10):
    filename = str(i) + filename_end
    platelayout180 = copy.copy(platelayout1)
    platelayout180.reverse() # I accidentally rotated the plate 180 degrees this day
    days[i-9] = prdl.import_platereader(filename, platelayout180)

for i in range(10,19):
    filename = str(i) + filename_end
    days[i-9] = prdl.import_platereader(filename,platelayout1)

os.chdir("../Experiment 1 Branch off July 18")

filename_end_branch = ' 2018 NS001 delta CAT mutation rate titrations.xlsx'

platelayout_branch = platelayout1

treatments = copy.deepcopy(platelayout_branch)
treatments.remove('blank1')
treatments.remove('blank2')
treatments.remove('blank3')
        
for i in range(0,8):
    filename = str(i+24) + ' August' + filename_end_branch
    days[i+10]=prdl.import_platereader(filename, platelayout_branch)

for i in range(8,17):
    filename = str(i-7) + ' September' + filename_end_branch
    days[i+10]=prdl.import_platereader(filename, platelayout_branch)

os.chdir("../Experiment 1 Branch 1 Continued 28 Sept")

filename_end_continue = ' 2018 NS001 delta CAT mutation rate titrations.xlsx'
        
for i in range(0,4):
    filename = str(i+27) + ' September' + filename_end_branch
    days[i+27]=prdl.import_platereader(filename, platelayout_branch)

for i in range(4,15):
    filename = str(i-3) + ' October' + filename_end_branch
    days[i+27] = prdl.import_platereader(filename, platelayout_branch)

The days are numbered 0 through 41. The replicates are named by relative mutation rate and replicate number. We call them High, HiMid, Mid, LoMid, and Low.

From lowest to highest mutation rate the corresponding inducer conditions of IPTG (micromolar) and aTc (ng/mL) are (2000,10), (2000,2), (2000,0), (50,0), and (0,0).

What do these curves typically look like? Let's see Day 1 and Day 41. We'll estimate the blank from the median of all blanks across all days. Earlier analysis showed a fairly consistent blank well from day to day, so to avoid overfitting, we'll use the same level of blank every day.

Growth curves on Day 1 vs Day 41

def estimate_blank(data):
    blanks = []
    for key in data:
        b1 = np.mean(data[key]['well_ODs']['blank1'])
        b2 = np.mean(data[key]['well_ODs']['blank2'])
        b3 = np.mean(data[key]['well_ODs']['blank3'])
        blanks.append(np.median(np.array([b1,b2,b3])))
    return np.median(np.array(blanks))

fig, axes = plt.subplots(2, 3, sharey='row', figsize=(18,12))
t = days[1]['times']
P = days[1]['well_ODs']['High 1'] - estimate_blank(days)
dP_dt = np.diff(P)/np.diff(t)
axes[0,0].plot(t, P);
axes[0,0].set_xlim([0,25]);
axes[0,0].set_xlabel('time (hours)', fontsize=18);
axes[0,0].set_ylabel('Optical Density (minus blank)', fontsize=18);
axes[0,0].text(.1, .9, 'A', transform=axes[0,0].transAxes, fontsize=18);
axes[1,0].plot(P[:-1], dP_dt, '*');
axes[1,0].set_xlim([-.03,.86]);
axes[1,0].set_xlabel('Optical Density (minus blank)', fontsize=18);
axes[1,0].set_ylabel('d(Optical Density)/dt', fontsize=18);
axes[1,0].text(.1, .9, 'D', transform=axes[1,0].transAxes, fontsize=18);

for sample in treatments:
    t = days[1]['times']
    P = days[1]['well_ODs'][sample] - estimate_blank(days)
    dP_dt = np.diff(P)/np.diff(t)
    axes[0,1].plot(t, P);
    axes[1,1].plot(P[:-1], dP_dt, '*');
axes[0,1].set_xlabel('time (hours)', fontsize=18);
axes[0,1].set_xlim([0,25]);
axes[0,1].text(.1, .9, 'B', transform=axes[0,1].transAxes, fontsize=18);
axes[1,1].set_xlim([-.03,.86]);
axes[1,1].set_xlabel('Optical Density (minus blank)', fontsize=18);
axes[1,1].text(.1, .9, 'E', transform=axes[1,1].transAxes, fontsize=18);

for sample in treatments:
    t = days[41]['times']
    P = days[41]['well_ODs'][sample] - estimate_blank(days)
    dP_dt = np.diff(P)/np.diff(t)
    axes[0,2].plot(t, P);
    axes[1,2].plot(P[:-1], dP_dt, '*');
axes[0,2].set_xlabel('time (hours)', fontsize=18);
axes[0,2].set_xlim([0,25]);
axes[0,2].text(.1, .9, 'C', transform=axes[0,2].transAxes, fontsize=18);
axes[1,2].set_xlim([-.03,.86]);
axes[1,2].set_xlabel('Optical Density (minus blank)', fontsize=18);
axes[1,2].text(.1, .9, 'F', transform=axes[1,2].transAxes, fontsize=18);

Extracting parameters describing the growth curves each day

There's clearly a great deal of regularity in the growth curves. An initial period where the signal to noise is too poor to get a good measurement, followed by exponential growth. Then a more complicated transition region. Then what looks like an exponential decay towards a saturating optical density.

So I'll analyze the data by extracting time to leave exponential, saturating optical density, and a curve fit for exponential phase

def exp(t, P0, r):
    return P0*np.exp(r*t)

def fit_exp(t, P, cut_below, cut_above):
    try:
        i1 = np.where(P<cut_below)[0][-1]+1
    except IndexError:
        i1 = 0
    i2 = np.where(P<cut_above)[0][-1]+1
    try:
        popt_exp, pcov_exp = spopt.curve_fit(exp, t[i1:i2], P[i1:i2], p0=(.002, 1))
    except RuntimeError:
        popt_exp = np.ones(3)*np.nan
        pcov_exp = np.ones((3,3))*np.nan
    return popt_exp, pcov_exp

def time_to_leave_exp(t, P, cut_above):
    i1 = np.where(P<cut_above)[0][-1]
    i2 = i1+1
    t1 = t[i1]
    t2 = t[i2]
    p1 = P[i1]
    p2 = P[i2]
    return t1 + (cut_above - p1)/(p2-p1)*(t2-t1)

def fit_over_experiments(data, days, samples, cut_below=.12, cut_above=.24):
    # set up the variables that will be returned with parameter estimates
    popt_exp = OrderedDict()
    pcov_exp = OrderedDict()
    K_est = OrderedDict()
    t_cut2 = OrderedDict()
    # estimate OD of empty well (blank) and adjust cutoffs by subtracting value of blank
    blank = estimate_blank(data)
    cut_below = cut_below - blank
    cut_above = cut_above - blank
    # estimate the parameters of exponential phase and the saturating OD
    for sample in samples:
        popt_exp[sample]=np.zeros((2,days))
        pcov_exp[sample]=np.zeros((2,2,days))
        K_est[sample]=np.zeros(days)
        t_cut2[sample]=np.zeros(days)
        for i in range(days):
            t = data[i]['times']
            P = data[i]['well_ODs'][sample] - blank
            pexp, covexp = fit_exp(t, P, cut_below, cut_above) 
            popt_exp[sample][:,i] = pexp
            pcov_exp[sample][:,:,i] = covexp
            K_est[sample][i] = np.mean(P[-10:])
            t_cut2[sample][i] = time_to_leave_exp(t, P, cut_above)
    return popt_exp, pcov_exp, K_est, t_cut2, blank

pexps, covexps, Ks, t_trans, blank = fit_over_experiments(days, 42, treatments)

What do these curve fits typically look like? A quick plot helps visualize how I'm extracting parameters.

fig, axes = plt.subplots(1, 2, figsize=(20,10))
dy = 32
smp = 'HiMid 1'
t = days[dy]['times']
P = days[dy]['well_ODs'][smp] - estimate_blank(days)
i1 = np.where(P<.12-estimate_blank(days))[0][-1]+1
i2 = np.where(P<.24-estimate_blank(days))[0][-1]+1
dP_dt = np.diff(P)/np.diff(t)
exp_phase = pexps[smp][0][dy]*np.exp(pexps[smp][1][dy]*t[i1:i2])
sat_OD = Ks[smp][dy]*np.ones(10)

axes[0].semilogy(t, P, 'b*');
axes[0].set_xlim([0,25]);
axes[0].set_xlabel('time (hours)', fontsize=18);
axes[0].set_ylabel('Optical Density (minus blank)', fontsize=18);
axes[0].text(.1, .9, 'A', transform=axes[0].transAxes, fontsize=18);
axes[0].axvline(t_trans[smp][dy],color='green');
axes[0].axhline(.12-estimate_blank(days), color='pink');
axes[0].axhline(.24-estimate_blank(days),color='red');
axes[0].semilogy(t[i1:i2], exp_phase, 'k');
axes[0].semilogy(t[-10:], sat_OD, 'y', linewidth=4);

axes[1].plot(P[:-10], dP_dt[:-9], '*');
axes[1].plot(P[-10:-1], dP_dt[-9:], 'y*');
axes[1].plot(P[i1:i2], pexps[smp][1][dy]*P[i1:i2],'k');
axes[1].set_xlim([-.03,.86]);
axes[1].set_xlabel('Optical Density (minus blank)', fontsize=18);
axes[1].set_ylabel('d(Optical Density)/dt', fontsize=18);
axes[1].text(.1, .9, 'B', transform=axes[1].transAxes, fontsize=18);
axes[1].axvline(.12-estimate_blank(days), color='pink');
axes[1].axvline(.24-estimate_blank(days),color='red');

Now I can plot the evolution of the exponential growth rate of the population over time, the evolution of the time to leave exponential, and the evolution of the saturating OD. Because there are 45 replicates, I show the data by mutation rate group. I plot slashes on the x-axis on the days where the experiment was interrupted by freezing and unfreezing. I also drop data from the first day, and from each day after an unfreezing because the initial population was not the same on these days as on other days and it affects the resulting parameters of the growth curve.

days_s = [i for i in range(1,42) if i != 27 and i !=10] # days that aren't the first day or after waking up from freezing
markers = ['o', 'v', '^', '<', '>', '+', 's', 'D', '*']
colors = ['C0', 'C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7', 'C8']

Evolution of the growth rate

[Return to Table of Contents](#Table-of-Contents)

It turns out that the growth rate does not evolve much over time. At least not in a way discernible from noise. Although some curves seem to show small decreases in their exponential phase growth over time.

plt.figure(figsize=(10,10))
param = 1
for i, sample in enumerate(['Low '+str(i) for i in range(1,10)]):
    plt.errorbar(np.arange(42)[days_s], pexps[sample][param][days_s], yerr=covexps[sample][param,param][days_s]**.5,
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of r at low mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('r', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.legend()
plt.ylim(.6,1.0);

plt.figure(figsize=(10,10))
param = 1
for i, sample in enumerate(['LoMid '+str(i) for i in range(1,10)]):
    plt.errorbar(np.arange(42)[days_s], pexps[sample][param][days_s], yerr=covexps[sample][param,param][days_s]**.5,
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of r at lomid mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('r', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.legend()
plt.ylim(.6,1.0);

plt.figure(figsize=(10,10))
param = 1
for i, sample in enumerate(['Mid '+str(i) for i in range(1,10)]):
    plt.errorbar(np.arange(42)[days_s], pexps[sample][param][days_s], yerr=covexps[sample][param,param][days_s]**.5,
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of r at mid mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('r', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.legend()
plt.ylim(.6,1.0);

plt.figure(figsize=(10,10))
param = 1
for i, sample in enumerate(['HiMid '+str(i) for i in range(1,10)]):
    plt.errorbar(np.arange(42)[days_s], pexps[sample][param][days_s], yerr=covexps[sample][param,param][days_s]**.5,
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of r at himid mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('r', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.legend()
plt.ylim(.6,1.0);

plt.figure(figsize=(10,10))
param = 1
for i, sample in enumerate(['High '+str(i) for i in range(1,10)]):
    plt.errorbar(np.arange(42)[days_s], pexps[sample][param][days_s], yerr=covexps[sample][param,param][days_s]**.5,
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of r at high mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('r', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.legend()
plt.ylim(.6,1.0);

We can graph the mean and coefficient of variation of the growth rate across all replicates at each mutation rate over time.

plt.figure(figsize=(10,10))
r_lows = np.array([pexps[sample][1,:] for sample in ['Low '+str(i) for i in range(1,10)]])
r_lowmids = np.array([pexps[sample][1,:] for sample in ['LoMid '+str(i) for i in range(1,10)]])
r_mids = np.array([pexps[sample][1,:] for sample in ['Mid '+str(i) for i in range(1,10)]])
r_highmids = np.array([pexps[sample][1,:] for sample in ['HiMid '+str(i) for i in range(1,10)]])
r_highs = np.array([pexps[sample][1,:] for sample in ['High '+str(i) for i in range(1,10)]])
plt.plot(np.arange(42)[days_s], np.mean(r_lows,axis=0)[days_s], linestyle='-', marker=markers[0], label='low mutation rate')
plt.plot(np.arange(42)[days_s], np.mean(r_lowmids,axis=0)[days_s], linestyle='-', marker=markers[2], label='lowmid mutation rate')
plt.plot(np.arange(42)[days_s], np.mean(r_mids,axis=0)[days_s], linestyle='-', marker=markers[6], label='mid mutation rate')
plt.plot(np.arange(42)[days_s], np.mean(r_highmids,axis=0)[days_s], linestyle='-', marker=markers[7], label='highmid mutation rate')
plt.plot(np.arange(42)[days_s], np.mean(r_highs,axis=0)[days_s], linestyle='-', marker=markers[8], label='high mutation rate')
plt.xlabel('day', fontsize=20);
plt.ylabel('mean growth rate', fontsize=20);
plt.ylim(.6,1.0);
plt.plot([26.5,27.5],[.58,.62],'k',clip_on=False)
plt.plot([26.5,27.5],[.98,1.02],'k',clip_on=False)
plt.plot([9.5,10.5],[.58,.62],'k',clip_on=False)
plt.plot([9.5,10.5],[.98,1.02],'k',clip_on=False)

[<matplotlib.lines.Line2D at 0x22669cd2c50>]

plt.figure(figsize=(10,10))
plt.plot(np.arange(42)[days_s], (np.std(r_lows,axis=0)/np.mean(r_lows,axis=0))[days_s], linestyle='-', marker=markers[0], label='low mutation rate')
plt.plot(np.arange(42)[days_s], (np.std(r_lowmids,axis=0)/np.mean(r_lowmids,axis=0))[days_s], linestyle='-', marker=markers[2], label='lowmid mutation rate')
plt.plot(np.arange(42)[days_s], (np.std(r_mids,axis=0)/np.mean(r_mids,axis=0))[days_s], linestyle='-', marker=markers[6], label='mid mutation rate')
plt.plot(np.arange(42)[days_s], (np.std(r_highmids,axis=0)/np.mean(r_highmids,axis=0))[days_s], linestyle='-', marker=markers[7], label='highmid mutation rate')
plt.plot(np.arange(42)[days_s], (np.std(r_highs,axis=0)/np.mean(r_highs,axis=0))[days_s], linestyle='-', marker=markers[8], label='high mutation rate')
plt.xlabel('day', fontsize=20);
plt.ylabel('coefficient of variation of growth rate', fontsize=20);
plt.ylim(0,.09);
plt.plot([26.5,27.5],[-.0025,.0025],'k',clip_on=False)
plt.plot([26.5,27.5],[.0875,.0925],'k',clip_on=False)
plt.plot([9.5,10.5],[-.0025,.0025],'k',clip_on=False)
plt.plot([9.5,10.5],[.0875,.0925],'k',clip_on=False)

[<matplotlib.lines.Line2D at 0x2266ae5d358>]

Everything looks pretty flat for r except for weird behavior right after unfreezing (expected) and around day 35 but only for low and lowmid mutation rates (not so expected). Let's graph some data from that day at those mutation rates to see if the problem is obvious.

plt.figure(figsize=(10,10))
all_days=list(range(36,37))
for i in all_days:
    for sample in ['Mid '+str(i) for i in range(1,10)]:
        t = days[i]['times']
        P = days[i]['well_ODs'][sample]
        plt.plot(t, P, '-', label = 'sample {}, day {}'.format(sample, i))
plt.ylabel('Optical Density');
plt.xlabel('Time (hrs.)');
plt.title('Typical curves over a day');
plt.legend()

<matplotlib.legend.Legend at 0x2266a551588>

The problem is obvious, and almost certainly with the measuring device, but I have no explanation for it. This is unfortunate.

Evolution of time to reach the transition point between exponential phase and after exponential phase

plt.figure(figsize=(10,10))
for i, sample in enumerate(['Low '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_s], t_trans[sample][days_s],
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of time to leave exponential at low mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('time to leave exponential phase (hours)', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.ylim((4.5,6.75))
plt.legend();

plt.figure(figsize=(10,10))
for i, sample in enumerate(['LoMid '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_s], t_trans[sample][days_s],
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of time to leave exponential at lowmid mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('time to leave exponential phase (hours)', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.ylim((4.5,6.75))
plt.legend();

plt.figure(figsize=(10,10))
for i, sample in enumerate(['Mid '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_s], t_trans[sample][days_s],
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of time to leave exponential at mid mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('time to leave exponential phase (hours)', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.ylim((4.5,6.75))
plt.legend();

plt.figure(figsize=(10,10))
for i, sample in enumerate(['HiMid '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_s], t_trans[sample][days_s],
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of time to leave exponential at himid mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('time to leave exponential phase (hours)', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.ylim((4.5,6.75))
plt.legend();

plt.figure(figsize=(10,10))
for i, sample in enumerate(['High '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_s], t_trans[sample][days_s],
                 linestyle='', marker=markers[i], color=colors[i], label=sample)
plt.title('evolution of time to leave exponential at high mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('time to leave exponential phase (hours)', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);
plt.ylim((4.5,6.75))
plt.legend();

plt.figure(figsize=(10,10))
trans_low = np.array([t_trans[sample] for sample in ['Low ' + str(i) for i in range(1,10)]])
trans_lomid = np.array([t_trans[sample] for sample in ['LoMid ' + str(i) for i in range(1,10)]])
trans_mid = np.array([t_trans[sample] for sample in ['Mid ' + str(i) for i in range(1,10)]])
trans_himid = np.array([t_trans[sample] for sample in ['HiMid ' + str(i) for i in range(1,10)]])
trans_hi = np.array([t_trans[sample] for sample in ['High ' + str(i) for i in range(1,10)]])
plt.plot(np.arange(0,42)[days_s], np.mean(trans_low,axis=0)[days_s], linestyle='-', marker=markers[0], label='low mutation rate')
plt.plot(np.arange(0,42)[days_s], np.mean(trans_lomid,axis=0)[days_s], linestyle='-', marker=markers[2], label='lowmid mutation rate')
plt.plot(np.arange(0,42)[days_s], np.mean(trans_mid,axis=0)[days_s], linestyle='-', marker=markers[6], label='mid mutation rate')
plt.plot(np.arange(0,42)[days_s], np.mean(trans_himid,axis=0)[days_s], linestyle='-', marker=markers[7], label='highmid mutation rate')
plt.plot(np.arange(0,42)[days_s], np.mean(trans_hi,axis=0)[days_s], linestyle='-', marker=markers[8], label='high mutation rate')
plt.ylim((4.65,6.55))
plt.plot([26.5,27.5],[4.6,4.7],'k',clip_on=False)
plt.plot([26.5,27.5],[6.5,6.6],'k',clip_on=False)
plt.plot([9.5,10.5],[4.6,4.7],'k',clip_on=False)
plt.plot([9.5,10.5],[6.5,6.6],'k',clip_on=False)
plt.xlabel('day', fontsize=20);
plt.ylabel('time to leave exponential phase (hours)', fontsize=20);

plt.figure(figsize=(10,10))
plt.plot(np.arange(0,42)[days_s], np.std(trans_low,axis=0)[days_s]/np.mean(trans_low,axis=0)[days_s],
         linestyle='-', marker=markers[0], label='low mutation rate')
plt.plot(np.arange(0,42)[days_s], np.std(trans_lomid,axis=0)[days_s]/np.mean(trans_lomid,axis=0)[days_s],
         linestyle='-', marker=markers[2], label='lowmid mutation rate')
plt.plot(np.arange(0,42)[days_s], np.std(trans_mid,axis=0)[days_s]/np.mean(trans_mid,axis=0)[days_s],
         linestyle='-', marker=markers[6], label='mid mutation rate')
plt.plot(np.arange(0,42)[days_s], np.std(trans_himid,axis=0)[days_s]/np.mean(trans_himid,axis=0)[days_s],
         linestyle='-', marker=markers[7], label='highmid mutation rate')
plt.plot(np.arange(0,42)[days_s], np.std(trans_hi,axis=0)[days_s]/np.mean(trans_hi,axis=0)[days_s],
         linestyle='-', marker=markers[8], label='high mutation rate')
plt.ylim(0,.11)
plt.plot([26.5,27.5],[-.005,.005],'k',clip_on=False)
plt.plot([26.5,27.5],[.105,.115],'k',clip_on=False)
plt.plot([9.5,10.5],[-.005,.005],'k',clip_on=False)
plt.plot([9.5,10.5],[.105,.115],'k',clip_on=False)
plt.xlabel('day', fontsize=20);
plt.ylabel('coefficient of variation of time to\nleave exponential phase (hours)', fontsize=20);

Evolution of saturating OD

[Return to Table of Contents](#Table-of-Contents)

The estimated saturating OD definitely evolves over time though at every set of mutation rates. There is a weird break in saturating OD between the first freeze and unfreeze though and I'm not sure what caused it. I'm pretty sure it is a systematic error and not a real effect. It looks like this.

plt.figure(figsize=(10,10))
for sample in treatments:
    plt.plot(np.arange(42)[days_s], Ks[sample][days_s], marker='s', linestyle='', label=sample)
plt.title('evolution of saturating OD at low mutation rates',fontsize=20);
plt.xlabel('day', fontsize=20);
plt.ylabel('saturating OD', fontsize=20);
plt.axvline(x=9.5);
plt.axvline(x=26.5);

I'll adjust for it by just dividing the values before the freeze by the ratio of their mean value compared to the mean value of two days after the freeze (skipping the wake up day of course).

fig = plt.figure(figsize=(12,12))
plt.ylim(.35,.95);
days_a = [i for i in range(11,42) if i != 27]
for i, sample in enumerate(['Low '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_a], Ks[sample][days_a],
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
days_b = list(range(1,10))
rough_mean_before_freeze = np.mean(Ks[sample][1:10])
for i, sample in enumerate(['Low '+str(i) for i in range(1,10)]):
    rough_mean_before_freeze = np.mean(Ks[sample][1:10])
    rough_mean_after_freeze = np.mean(Ks[sample][11:13])
    plt.plot(np.arange(42)[days_b], Ks[sample][days_b]*rough_mean_after_freeze/rough_mean_before_freeze,
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
plt.plot([26.5,27.5],[.33,.37],'k',clip_on=False)
plt.plot([26.5,27.5],[.93,.97],'k',clip_on=False)
plt.plot([9.5,10.5],[.33,.37],'k',clip_on=False)
plt.plot([9.5,10.5],[.93,.97],'k',clip_on=False)
plt.xlabel('day', fontsize=20);
plt.ylabel('saturating OD', fontsize=20);

At low mutation rates, we see the saturating OD doesn't evolve much until about midway, at which replicates at first tend to have increases in carrying capacity. However after a bit of time, carrying capacity suddenly plummets. We'll see a similar pattern at higher mutation rates except it occurs quicker.

fig = plt.figure(figsize=(12,12))
plt.ylim(.35,.95);
for i, sample in enumerate(['LoMid '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_a], Ks[sample][days_a],
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
rough_mean_before_freeze = np.mean(Ks[sample][1:10])
for i, sample in enumerate(['LoMid '+str(i) for i in range(1,10)]):
    rough_mean_before_freeze = np.mean(Ks[sample][1:10])
    rough_mean_after_freeze = np.mean(Ks[sample][11:13])
    plt.plot(np.arange(42)[days_b], Ks[sample][days_b]*rough_mean_after_freeze/rough_mean_before_freeze,
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
plt.plot([26.5,27.5],[.33,.37],'k',clip_on=False)
plt.plot([26.5,27.5],[.93,.97],'k',clip_on=False)
plt.plot([9.5,10.5],[.33,.37],'k',clip_on=False)
plt.plot([9.5,10.5],[.93,.97],'k',clip_on=False)
plt.xlabel('day', fontsize=20);
plt.ylabel('saturating OD', fontsize=20);

fig = plt.figure(figsize=(12,12))
plt.ylim(.35,.95);
for i, sample in enumerate(['Mid '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_a], Ks[sample][days_a],
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
rough_mean_before_freeze = np.mean(Ks[sample][1:10])
for i, sample in enumerate(['Mid '+str(i) for i in range(1,10)]):
    rough_mean_before_freeze = np.mean(Ks[sample][1:10])
    rough_mean_after_freeze = np.mean(Ks[sample][11:13])
    plt.plot(np.arange(42)[days_b], Ks[sample][days_b]*rough_mean_after_freeze/rough_mean_before_freeze,
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
plt.plot([26.5,27.5],[.33,.37],'k',clip_on=False)
plt.plot([26.5,27.5],[.93,.97],'k',clip_on=False)
plt.plot([9.5,10.5],[.33,.37],'k',clip_on=False)
plt.plot([9.5,10.5],[.93,.97],'k',clip_on=False)
plt.xlabel('day', fontsize=20);
plt.ylabel('saturating OD', fontsize=20);

fig = plt.figure(figsize=(12,12))
plt.ylim(.35,.95);
for i, sample in enumerate(['HiMid '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_a], Ks[sample][days_a],
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
rough_mean_before_freeze = np.mean(Ks[sample][1:10])
for i, sample in enumerate(['HiMid '+str(i) for i in range(1,10)]):
    rough_mean_before_freeze = np.mean(Ks[sample][1:10])
    rough_mean_after_freeze = np.mean(Ks[sample][11:13])
    plt.plot(np.arange(42)[days_b], Ks[sample][days_b]*rough_mean_after_freeze/rough_mean_before_freeze,
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
plt.plot([26.5,27.5],[.33,.37],'k',clip_on=False)
plt.plot([26.5,27.5],[.93,.97],'k',clip_on=False)
plt.plot([9.5,10.5],[.33,.37],'k',clip_on=False)
plt.plot([9.5,10.5],[.93,.97],'k',clip_on=False)
plt.xlabel('day', fontsize=20);
plt.ylabel('saturating OD', fontsize=20);

fig = plt.figure(figsize=(12,12))
plt.ylim(.35,.95);
for i, sample in enumerate(['High '+str(i) for i in range(1,10)]):
    plt.plot(np.arange(42)[days_a], Ks[sample][days_a],
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
rough_mean_before_freeze = np.mean(Ks[sample][1:10])
for i, sample in enumerate(['High '+str(i) for i in range(1,10)]):
    rough_mean_before_freeze = np.mean(Ks[sample][1:10])
    rough_mean_after_freeze = np.mean(Ks[sample][11:13])
    plt.plot(np.arange(42)[days_b], Ks[sample][days_b]*rough_mean_after_freeze/rough_mean_before_freeze,
                 linestyle='', marker=markers[i], markersize=7, color=colors[i], label=sample)
plt.plot([26.5,27.5],[.33,.37],'k',clip_on=False)
plt.plot([26.5,27.5],[.93,.97],'k',clip_on=False)
plt.plot([9.5,10.5],[.33,.37],'k',clip_on=False)
plt.plot([9.5,10.5],[.93,.97],'k',clip_on=False)
plt.xlabel('day', fontsize=20);
plt.ylabel('saturating OD', fontsize=20);

At high mutation rates, carrying capacity decreases much faster. We can see this side by side if we average across replicates at a mutation rate.

K_lows = np.array([Ks[sample] for sample in ['Low '+str(i) for i in range(1,10)]])
K_lowmids = np.array([Ks[sample] for sample in ['LoMid '+str(i) for i in range(1,10)]])
K_mids = np.array([Ks[sample] for sample in ['Mid '+str(i) for i in range(1,10)]])
K_highmids = np.array([Ks[sample] for sample in ['HiMid '+str(i) for i in range(1,10)]])
K_highs = np.array([Ks[sample] for sample in ['High '+str(i) for i in range(1,10)]])

K_lows_adj = np.copy(K_lows)
K_lows_adj[:,1:10]=K_lows_adj[:,1:10]/np.mean(K_lows_adj[:,1:10],axis=1)*np.mean(K_lows_adj[:,11:13],axis=1)
K_lowmids_adj = np.copy(K_lowmids)
K_lowmids_adj[:,1:10]=K_lowmids_adj[:,1:10]/np.mean(K_lowmids_adj[:,1:10],axis=1)*np.mean(K_lowmids_adj[:,11:13],axis=1)
K_mids_adj = np.copy(K_mids)
K_mids_adj[:,1:10]=K_mids_adj[:,1:10]/np.mean(K_mids_adj[:,1:10],axis=1)*np.mean(K_mids_adj[:,11:13],axis=1)
K_highmids_adj = np.copy(K_highmids)
K_highmids_adj[:,1:10]=K_highmids_adj[:,1:10]/np.mean(K_highmids_adj[:,1:10],axis=1)*np.mean(K_highmids_adj[:,11:13],axis=1)
K_highs_adj = np.copy(K_highs)
K_highs_adj[:,1:10]=K_highs_adj[:,1:10]/np.mean(K_highs_adj[:,1:10],axis=1)*np.mean(K_highs_adj[:,11:13],axis=1)

plt.figure(figsize=(10,10))
plt.ylim(.35,.95)
plt.plot(np.arange(42)[days_s], np.mean(K_lows_adj,axis=0)[days_s],
         linestyle='-', marker=markers[0], label='low mutation rate')
plt.plot(np.arange(42)[days_s], np.mean(K_lowmids_adj,axis=0)[days_s],
         linestyle='-', marker=markers[2], label='lowmid mutation rate')
plt.plot(np.arange(42)[days_s], np.mean(K_mids_adj,axis=0)[days_s],
         linestyle='-', marker=markers[6], label='mid mutation rate')
plt.plot(np.arange(42)[days_s], np.mean(K_highmids_adj,axis=0)[days_s],
         linestyle='-', marker=markers[7], label='highmid mutation rate')
plt.plot(np.arange(42)[days_s], np.mean(K_highs_adj,axis=0)[days_s],
         linestyle='-', marker=markers[8], label='high mutation rate')

plt.plot([26.5,27.5],[.33,.37],'k',clip_on=False)
plt.plot([26.5,27.5],[.93,.97],'k',clip_on=False)
plt.plot([9.5,10.5],[.33,.37],'k',clip_on=False)
plt.plot([9.5,10.5],[.93,.97],'k',clip_on=False)
plt.xlabel('day', fontsize=20);
plt.ylabel('mean saturating OD', fontsize=20);

Another important thing to note is that the coefficient of variation of saturating OD for a given mutation rate condition increases over time.

plt.figure(figsize=(10,10))
plt.plot(np.arange(42)[days_s], (np.std(K_lows_adj,axis=0)/np.mean(K_lows_adj,axis=0))[days_s],
         linestyle='-', marker=markers[0], label='low mutation rate')
plt.plot(np.arange(42)[days_s], (np.std(K_lowmids_adj,axis=0)/np.mean(K_lowmids_adj,axis=0))[days_s],
         linestyle='-', marker=markers[2], label='lowmid mutation rate')
plt.plot(np.arange(42)[days_s], (np.std(K_mids_adj,axis=0)/np.mean(K_mids_adj,axis=0))[days_s],
         linestyle='-', marker=markers[6], label='mid mutation rate')
plt.plot(np.arange(42)[days_s], (np.std(K_highmids_adj,axis=0)/np.mean(K_highmids_adj,axis=0))[days_s],
         linestyle='-', marker=markers[7], label='highmid mutation rate')
plt.plot(np.arange(42)[days_s], (np.std(K_highs_adj,axis=0)/np.mean(K_highs_adj,axis=0))[days_s],
         linestyle='-', marker=markers[8], label='high mutation rate')
plt.xlabel('day', fontsize=20);
plt.ylabel('coefficient of variation of saturating OD', fontsize=20);
plt.ylim(0,.2)
plt.plot([26.5,27.5],[-.005,.005],'k',clip_on=False)
plt.plot([26.5,27.5],[.195,.205],'k',clip_on=False)
plt.plot([9.5,10.5],[-.005,.005],'k',clip_on=False)
plt.plot([9.5,10.5],[.195,.205],'k',clip_on=False)

[<matplotlib.lines.Line2D at 0x2266b430240>]

Blanks and cross-contamination

Blanks are fairly uniform across the experiment. We saw a few instances of cross-contamination of blanks in the first six days of the experiment. However, after using tape to cover wells during transfer and moving as carefully as possible when pipetting, we ceased to see cross contamination for the remaining 35 days. Thus we think cross-contamination did not unduly affect the results.

plt.figure(figsize=(10,10))
for dy in range(42):
    plt.plot(days[dy]['times'],days[dy]['well_ODs']['blank1'],label='blank1, day' + str(dy))
    plt.plot(days[dy]['times'],days[dy]['well_ODs']['blank2'],label='blank2, day' + str(dy))
    plt.plot(days[dy]['times'],days[dy]['well_ODs']['blank3'],label='blank3, day' + str(dy))

plt.figure(figsize=(10,10))
for dy in range(2,6):
    plt.plot(days[dy]['times'],days[dy]['well_ODs']['blank1'],label='blank1, day' + str(dy))
    plt.plot(days[dy]['times'],days[dy]['well_ODs']['blank2'],label='blank2, day' + str(dy))
    plt.plot(days[dy]['times'],days[dy]['well_ODs']['blank3'],label='blank3, day' + str(dy))
plt.legend()

<matplotlib.legend.Legend at 0x2266a05f940>

fig, axes = plt.subplots(1, 2, figsize=(20,10))
dy = 2
smp = 'blank1'
t = days[dy]['times']
P = days[dy]['well_ODs'][smp] - estimate_blank(days)
i1 = np.where(P<.12-estimate_blank(days))[0][-1]+1
i2 = np.where(P<.24-estimate_blank(days))[0][-1]+1
dP_dt = np.diff(P)/np.diff(t)

axes[0].semilogy(t, P, 'b*');
axes[0].set_xlim([0,25]);
axes[0].set_xlabel('time (hours)', fontsize=18);
axes[0].set_ylabel('Optical Density (minus blank)', fontsize=18);
axes[0].text(.1, .9, 'A', transform=axes[0].transAxes, fontsize=18);
axes[0].axhline(.12-estimate_blank(days), color='pink');
axes[0].axhline(.24-estimate_blank(days),color='red');
axes[0].axvline(t[i2],color='green');

axes[1].plot(P[:-10], dP_dt[:-9], '*');
axes[1].plot(P[-10:-1], dP_dt[-9:], 'y*');
axes[1].set_xlim([-.03,.86]);
axes[1].set_xlabel('Optical Density (minus blank)', fontsize=18);
axes[1].set_ylabel('d(Optical Density)/dt', fontsize=18);
axes[1].text(.1, .9, 'B', transform=axes[1].transAxes, fontsize=18);
axes[1].axvline(.12-estimate_blank(days), color='pink');
axes[1].axvline(.24-estimate_blank(days),color='red');

There is a small spike in the early OD of many blanks that is part of what ruins our chances of getting a good measure of lag phase. It's easiest to see by averaging across blank1 for the latter half of the experiment. But it doesn't always happen, which makes correcting for it unwise. Better to do experiment with a higher starting OD so lag phase can be measured.

blnk_arr = np.array([days[dy]['well_ODs'][smp] for dy in range(15,42) for smp in ['blank1']
                     if len(days[dy]['times'])==128])
print(blnk_arr.shape)

(12, 128)

plt.plot(days[40]['times'][:],np.mean(blnk_arr,axis=0)[:],'*')

[<matplotlib.lines.Line2D at 0x2266b45a7f0>]