Analysis of growth data in the efflux experiment

(c) 2024 Justin Bois. This work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.

This document was prepared at Caltech with support financial support from the Donna and Benjamin M. Rosen Bioengineering Center.

This tutorial was generated from an Jupyter notebook. You can download the notebook here.

---

import numpy as np
import pandas as pd
import scipy.optimize

import iqplot

import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook()

BokehJS 3.4.1 successfully loaded.

---

In this tutorial, we will explore how to analyze the data from the efflux pump experiment. We will then develop some tools for parsing the data set and obtaining growth rates from the curves. Then, we will learn how to compute the growth rate for a given set of parameters and antibiotic concentration. This is necessary to be able to perform a regression to obtain an estimate the parameter $q_{\mathrm{ex}, 0}^\mathrm{max}$ for each strain.

Data wrangling

The data will again come from our plate reader in the same format as for the E. coli growth experiment. You can use the parser from that module to load in the data into a data frame. Remember that you need to specify the filename, header row, and last row of the data to do so. Here is that parser for convenience.

def read_plate_reader(fname, header_row, last_row):
    """Read in time series data from Bi 1x plate reader.

    Parameters
    ----------
    fname : str
        Path to file outputted from plate reader.
    header_row : int
        Number of row containing header column (zero-indexed)
    last_row : int
        Last row of the file containing the time series data (zero-indexed)

    Returns
    -------
    output : DataFrame
        Pandas DataFrame containing time series data.
    """
    # Find out how many rows in input file
    with open(fname, "rb") as fp:
        n_rows = len(fp.readlines())
    
    # How many lines in footer to skip
    skipfooter = n_rows - last_row - 1

    # Read in data frame
    df = pd.read_csv(
        fname,
        sep="\t",
        skiprows=header_row,
        skipfooter=skipfooter,
        engine="python",
        encoding='ISO-8859-1',
    )

    # Rename temperature column
    df = df.rename(columns={df.columns[df.columns.str.contains('T°')][0]: 'temperature (deg C)'})

    # Parse the time column
    df['time (s)'] = pd.to_timedelta(df['Time']).dt.total_seconds()

    # Set start time to zero
    df['time (s)'] -= df['time (s)'].min()

    # Time in units of hours
    df['time (hr)'] = df['time (s)'] / 3600

    return df

Well metadata

In this experiment, we will have two data sets to load in, one from Group A, and one from Group B. You will perform regressions to get growth rates on each data set separately. That is, you will determine $\lambda_0$ and $\lambda$ for each time course, reporting relative growth rate and $\lambda/\lambda_0$, where the $\lambda_0$ is determined from the data set from the same plate. That is, you should have a $\lambda_0$ for Group A and a separate $\lambda_0$ for Group B.

For this demonstration, we will work with one data set from a dry run from this experiment in the file 20240527_tet_test.txt. The analysis is analogous for a second data set. As in the E. coli growth module, we need to specify the header row and the last row of data.

# This is a dry run of this experiment
fname_A = '20240527_tet_test.txt'
header_row_A = 53
last_row_A = 223

# Load in the data frame for "Group A" (which is just the dry run in this case)
dfA = read_plate_reader(fname_A, header_row_A, last_row_A)

It is also useful to have a dictionary with the metadata for each well in the plates, in this case the strain and tetracycline concentration. We have two dictionaries, one for the plate from Group A and one from Group B. Below are dictionaries for the metadata for this experiment. The abbreviated strain is followed by the tetracycline concentration in µM, separated by a dash. Blank wells are labeled as "blank."

well_metadata_A = dict(
    A1="MG-0.0",
    A2="WT-0.21",
    A3="MG-0.0",
    A4="WT-0.90",
    A5="MG-0.0",
    A6="WT-3.90",
    A7="MG-0.0",
    A8="WT-16.90",
    A9="MG-0.0",
    A10="WT-73.15",
    A11="MG-0.0",
    A12="WT-208.33",
    B1="3.19-0.0",
    B2="UV5-0.21",
    B3="3.19-0.0",
    B4="UV5-0.90",
    B5="3.19-0.0",
    B6="UV5-3.90",
    B7="3.19-0.0",
    B8="UV5-16.90",
    B9="3.19-0.0",
    B10="UV5-73.15",
    B11="3.19-0.0",
    B12="UV5-208.33",
    C1="IW-0.0",
    C2="MG-0.21",
    C3="IW-0.0",
    C4="MG-0.90",
    C5="IW-0.0",
    C6="MG-3.90",
    C7="IW-0.0",
    C8="MG-16.90",
    C9="IW-0.0",
    C10="MG-73.15",
    C11="IW-0.0",
    C12="MG-208.33",
    D1="WT-0.0",
    D2="3.19-0.21",
    D3="WT-0.0",
    D4="3.19-0.90",
    D5="WT-0.0",
    D6="3.19-3.90",
    D7="WT-0.0",
    D8="3.19-16.90",
    D9="WT-0.0",
    D10="3.19-73.15",
    D11="WT-0.0",
    D12="3.19-208.33",
    E1="UV5-0.0",
    E2="IW-0.21",
    E3="UV5-0.0",
    E4="IW-0.90",
    E5="UV5-0.0",
    E6="IW-3.90",
    E7="UV5-0.0",
    E8="IW-16.90",
    E9="UV5-0.0",
    E10="IW-73.15",
    E11="UV5-0.0",
    E12="IW-208.33",
    F1="MG-0.21",
    F2="WT-0.21",
    F3="MG-0.90",
    F4="WT-0.90",
    F5="MG-3.90",
    F6="WT-3.90",
    F7="MG-16.90",
    F8="WT-16.90",
    F9="MG-73.15",
    F10="WT-73.15",
    F11="MG-208.33",
    F12="WT-208.33",
    G1="3.19-0.21",
    G2="UV5-0.21",
    G3="3.19-0.90",
    G4="UV5-0.90",
    G5="3.19-3.90",
    G6="UV5-3.90",
    G7="3.19-16.90",
    G8="UV5-16.90",
    G9="3.19-73.15",
    G10="UV5-73.15",
    G11="3.19-208.33",
    G12="UV5-208.33",
    H1="IW-0.21",
    H2="blank",
    H3="IW-0.90",
    H4="blank",
    H5="IW-3.90",
    H6="blank",
    H7="IW-16.90",
    H8="blank",
    H9="IW-73.15",
    H10="blank",
    H11="IW-208.33",
    H12="blank",
)

well_metadata_B = dict(
    A1="MG-0.0",
    A2="WT-0.44",
    A3="MG-0.0",
    A4="WT-1.88",
    A5="MG-0.0",
    A6="WT-8.13",
    A7="MG-0.0",
    A8="WT-35.17",
    A9="MG-0.0",
    A10="WT-152.19",
    A11="MG-0.0",
    A12="WT-312.50",
    B1="3.19-0.0",
    B2="UV5-0.44",
    B3="3.19-0.0",
    B4="UV5-1.88",
    B5="3.19-0.0",
    B6="UV5-8.13",
    B7="3.19-0.0",
    B8="UV5-35.17",
    B9="3.19-0.0",
    B10="UV5-152.19",
    B11="3.19-0.0",
    B12="UV5-312.50",
    C1="IW-0.0",
    C2="MG-0.44",
    C3="IW-0.0",
    C4="MG-1.88",
    C5="IW-0.0",
    C6="MG-8.13",
    C7="IW-0.0",
    C8="MG-35.17",
    C9="IW-0.0",
    C10="MG-152.19",
    C11="IW-0.0",
    C12="MG-312.50",
    D1="WT-0.0",
    D2="3.19-0.44",
    D3="WT-0.0",
    D4="3.19-1.88",
    D5="WT-0.0",
    D6="3.19-8.13",
    D7="WT-0.0",
    D8="3.19-35.17",
    D9="WT-0.0",
    D10="3.19-152.19",
    D11="WT-0.0",
    D12="3.19-312.50",
    E1="UV5-0.0",
    E2="IW-0.44",
    E3="UV5-0.0",
    E4="IW-1.88",
    E5="UV5-0.0",
    E6="IW-8.13",
    E7="UV5-0.0",
    E8="IW-35.17",
    E9="UV5-0.0",
    E10="IW-152.19",
    E11="UV5-0.0",
    E12="IW-312.50",
    F1="MG-0.44",
    F2="WT-0.44",
    F3="MG-1.88",
    F4="WT-1.88",
    F5="MG-8.13",
    F6="WT-8.13",
    F7="MG-35.17",
    F8="WT-35.17",
    F9="MG-152.19",
    F10="WT-152.19",
    F11="MG-312.50",
    F12="WT-312.50",
    G1="3.19-0.44",
    G2="UV5-0.44",
    G3="3.19-1.88",
    G4="UV5-1.88",
    G5="3.19-8.13",
    G6="UV5-8.13",
    G7="3.19-35.17",
    G8="UV5-35.17",
    G9="3.19-152.19",
    G10="UV5-152.19",
    G11="3.19-312.50",
    G12="UV5-312.50",
    H1="IW-0.44",
    H2="blank",
    H3="IW-1.88",
    H4="blank",
    H5="IW-8.13",
    H6="blank",
    H7="IW-35.17",
    H8="blank",
    H9="IW-152.19",
    H10="blank",
    H11="IW-312.50",
    H12="blank",
)

Tidy data

To more easily work with this data set, we will put it in tidy format, where each row in the data frame is a single absorbance measurement and each column is a property associated with that measurement, such as time, time stamp, temperature, well, well metadata, and, of course absorbance. The function below converts the data to tidy format given the well metadata.

def conc(well, well_metadata):
    """Return abx concentration for a given well"""
    strain_conc = well_metadata[well]
    if strain_conc == "blank":
        return 0.0
    else:
        ind = strain_conc.find("-")
        return float(strain_conc[ind + 1 :])


def strain(well, well_metadata):
    """Return strain for a given well"""
    strain_conc = well_metadata[well]
    if strain_conc == "blank":
        return "blank"
    else:
        ind = strain_conc.find("-")
        return strain_conc[:ind]


def tidy(df, well_metadata, group):
    """Convert data frame to tidy format."""
    # Convert to tall format
    df_tidy = df.melt(
        id_vars=["Time", "temperature (deg C)", "time (s)", "time (hr)"],
        var_name="well",
        value_name="absorbance",
    )

    # Add metadata columns
    df_tidy["conc"] = df_tidy["well"].apply(conc, args=(well_metadata,))
    df_tidy["strain"] = df_tidy["well"].apply(strain, args=(well_metadata,))
    df_tidy["group"] = group

    df_tidy = df_tidy.sort_values(
        by=["group", "strain", "conc", "well", "time (hr)"], ignore_index=True
    )

    return df_tidy

We can put it to use to tidy our data frame.

dfA = tidy(dfA, well_metadata_A, "A")

# Take a look
dfA.head()

	Time	temperature (deg C)	time (s)	time (hr)	well	absorbance	conc	strain	group
0	0:05:10	37.0	0.0	0.0	B1	0.085	0.0	3.19	A
1	0:11:10	37.1	360.0	0.1	B1	0.085	0.0	3.19	A
2	0:17:10	37.1	720.0	0.2	B1	0.085	0.0	3.19	A
3	0:23:10	37.0	1080.0	0.3	B1	0.085	0.0	3.19	A
4	0:29:10	37.0	1440.0	0.4	B1	0.085	0.0	3.19	A

For simplicity, I have labeled the antibiotic concentration column "conc", but recall that the concentration is in units of µM.

Checking the blanks

As before, we can check the blanks. To do this, we first pull out the part of the data frame that contains blanks. We then use Pandas's convenient groupby operation to group the data frame by well. We can iterate over the GroupBy object with a for loop, which gives us a group name (which will be the well in this case) and a sub-data frame that has just the data for that well. We can then add the trace of that blank to the plot.

p = bokeh.plotting.figure(
    frame_width=350,
    frame_height=200,
    x_axis_label='time (hr)',
    y_axis_label='absorbance'
)

# Pull out only entries for blanks
dfA_blank = dfA.loc[dfA['strain']=='blank', :]

# Groupby and iterate, adding to plot
for group, sub_df in dfA_blank.groupby('well'):
    p.line(sub_df['time (hr)'], sub_df['absorbance'])

bokeh.io.show(p)

The blanks are all really consistent. We can therefore define our blank concentration to be the average of all blank measurement.

blank_absorbance = dfA_blank['absorbance'].mean()

We can now add a column to the data frame for the OD600, which is the background substracted absorbance.

dfA['OD600'] = dfA['absorbance'] - blank_absorbance

Make plots

We can now use this kind of groupby structure to make all kinds of plots. For example, we can plot growth curves for each strain without antibiotic.

p = bokeh.plotting.figure(
    frame_width=350,
    frame_height=200,
    x_axis_label="time (hr)",
    y_axis_label="OD600",
    y_axis_type="log",
    toolbar_location="above",
)

# Color scheme
strains = ["MG", "3.19", "IW", "WT", "UV5"]
colors = {strain: color for strain, color in zip(strains, bokeh.palettes.Category10_5)}

# Pull out only entries for zero antibiotic concentration (and not blanks)
inds = (dfA["conc"] == 0.0) & (dfA["strain"] != "blank")
dfA_0abx = dfA.loc[inds, :]

# Items we will place in the legend
legend_items = []

# Groupby and iterate, adding to plot
for strain, sub_df in dfA_0abx.groupby("strain"):
    # Lines for a given strain
    lines = []

    # Get each well for the strain and add a line to the plot
    for well, subsub_df in sub_df.groupby("well"):
        lines.append(
            p.line(subsub_df["time (hr)"], subsub_df["OD600"], color=colors[strain])
        )

    # Add the lines to the legend
    legend_items.append((strain, lines))

# Create the legend from the items
legend = bokeh.models.Legend(items=legend_items, click_policy="hide")

# Add the legend to the plot
p.add_layout(legend, "right")

bokeh.io.show(p)

We could also plot the growth for all concentrations of tetracycline for the MG1655 strain.

p = bokeh.plotting.figure(
    frame_width=350,
    frame_height=200,
    x_axis_label="time (hr)",
    y_axis_label="OD600",
    # y_axis_type="log",
    toolbar_location="above",
)

# Color scheme
concs = np.sort(dfA['conc'].unique())
colors = {conc: color for conc, color in zip(concs, bokeh.palettes.Blues9)}

# Pull out only entries for MG1655
dfA_MG = dfA.loc[dfA['strain'] == 'MG', :]

# Items we will place in the legend
legend_items = []

# Groupby and iterate, adding to plot
for conc, sub_df in dfA_MG.groupby("conc"):
    # Lines for a given strain
    lines = []

    # Get each well for the strain and add a line to the plot
    for well, subsub_df in sub_df.groupby("well"):
        lines.append(
            p.line(subsub_df["time (hr)"], subsub_df["OD600"], color=colors[conc])
        )

    # Add the lines to the legend
    legend_items.append((str(conc), lines))

# Create the legend from the items
legend = bokeh.models.Legend(items=legend_items, click_policy="hide")

# Add the legend to the plot
p.add_layout(legend, "right")

bokeh.io.show(p)

Apparently, above an antibiotic concentration of 0.9, the cells do not grow.

Obtaining growth rates

We can use the groupby operations to perform curve fits as well. We need to modify our fit_growth() function from the E. coli growth module for use here to simply return the growth rate (we do not care about the initial absorbance).

def exp_growth(t, A0, lam):
    return A0 * np.exp(lam * t)


def fit_growth(df, t_start, t_end):
    """Obtain estimates for A0 and λ for a given sub-trace"""
    # If no time window, return zero
    if t_start == t_end:
        return 0.0
        
    # Parameter guesses
    A0_guess = 0.01
    lam_guess = 1.0  # 1/hr

    # Pull out data
    sub_df = df.loc[(df["time (hr)"] >= t_start) & (df["time (hr)"] <= t_end), :]
    t = sub_df["time (hr)"].values
    A = sub_df['OD600'].values

    # Perform optimization
    popt, _ = scipy.optimize.curve_fit(exp_growth, t, A, p0=(A0_guess, lam_guess))

    # Just return growth rate 
    return popt[1]

Now that we have our curve fitting function in hand, we can, for example, get growth rates for each strain in the absence of antibiotic (an important calculation; these are the $\lambda_0$ values). It turns out that Pandas GroupBy objects are really slick for doing this automatically using the apply() method. We pass in the function we want to call, and then specify other arguments to the function.

Looking at the plot, the exponential growth phase seems to be between 6 and 10 hours.

# Pull out only entries for zero antibiotic concentration (and not blanks)
inds = (dfA["conc"] == 0.0) & (dfA["strain"] != "blank")
dfA_0abx = dfA.loc[inds, :]

# Include include_groups=False to silence deprecation warning
# reset_index() to convert to dataframe
growth_rates_0abx = (
    dfA_0abx.groupby(["strain", "well"])
    .apply(fit_growth, t_start=6, t_end=10, include_groups=False)
    .reset_index(name='growth rate (1/hr)')
)

# Take a look
growth_rates_0abx

	strain	well	growth rate (1/hr)
0	3.19	B1	0.506579
1	3.19	B11	0.439572
2	3.19	B3	0.565789
3	3.19	B5	0.545540
4	3.19	B7	0.574230
5	3.19	B9	0.586322
6	IW	C1	0.531716
7	IW	C11	0.447853
8	IW	C3	0.533252
9	IW	C5	0.537041
10	IW	C7	0.535596
11	IW	C9	0.479817
12	MG	A1	0.512562
13	MG	A11	0.515497
14	MG	A3	0.487880
15	MG	A5	0.497014
16	MG	A7	0.471504
17	MG	A9	0.491036
18	UV5	E1	0.527472
19	UV5	E11	0.516805
20	UV5	E3	0.514496
21	UV5	E5	0.502841
22	UV5	E7	0.499688
23	UV5	E9	0.463137
24	WT	D1	0.509833
25	WT	D11	0.458630
26	WT	D3	0.499100
27	WT	D5	0.479537
28	WT	D7	0.488544
29	WT	D9	0.492815

We now have growth rates for each curve without antibiotics! We can make a plot of these to make sure there are no outliers.

bokeh.io.show(iqplot.strip(growth_rates_0abx, q='growth rate (1/hr)', cats='strain'))

The IW and 3.19 strains both have two growth rates that are low, but they are only off by about 20%. We can take $\lambda_0$ to be the median of our measurements (the median being immune to outliers) to get our $\lambda_0$ values.

lam0 = (
    growth_rates_0abx.groupby("strain")["growth rate (1/hr)"]
    .apply(np.median)
)

# Take a look
lam0

strain
3.19    0.555665
IW      0.532484
MG      0.494025
UV5     0.508669
WT      0.490680
Name: growth rate (1/hr), dtype: float64

We can similarly obtain growth rates for all concentrations of antibiotics in the MG1655 strain. We need to refer to the plot to give time windows. For tetracycline concentration above 1 µM, we saw no growth, so we do not need to fit those, and we can put [0, 0] for the time window. We have already updated the fitting function to return zero for the growth rate if there is a time window of zero width.

# Pull out only entries MG1655 with nonzero ABX concentration
dfA_MG = dfA.loc[(dfA["strain"] == "MG") & (dfA['conc'] > 0), :]

# Time windows based on plots
time_windows = {
    0.21: [6, 10.5],
    0.9: [14, 17],
    3.9: [0, 0],
    16.9: [0, 0],
    73.15: [0, 0],
    208.33: [0, 0],
}

# Compute curve fits (have to loop this time because we need time windows)
# Instantiate dictionary with results
growth_rates_MG = {'conc': [], 'well': [], 'growth rate (1/hr)': []}

# Loop through groupby object to do curve fits
for (conc, well), sub_df in dfA_MG.groupby(['conc', 'well']):
    # Update concentration and well
    growth_rates_MG['conc'].append(conc)
    growth_rates_MG['well'].append(well)

    # Obtain growth rate
    lam = fit_growth(sub_df, t_start=time_windows[conc][0], t_end=time_windows[conc][1])
    growth_rates_MG['growth rate (1/hr)'].append(lam)

# Convert to data frame
growth_rates_MG = pd.DataFrame(growth_rates_MG)

# Take a look
growth_rates_MG

	conc	well	growth rate (1/hr)
0	0.21	C2	0.383576
1	0.21	F1	0.428159
2	0.90	C4	0.242548
3	0.90	F3	0.180286
4	3.90	C6	0.000000
5	3.90	F5	0.000000
6	16.90	C8	0.000000
7	16.90	F7	0.000000
8	73.15	C10	0.000000
9	73.15	F9	0.000000
10	208.33	C12	0.000000
11	208.33	F11	0.000000

We can then compute the growth rate as the average.

lam_MG = growth_rates_MG.groupby('conc')['growth rate (1/hr)'].mean().reset_index()

# Take a look
lam_MG

	conc	growth rate (1/hr)
0	0.21	0.405868
1	0.90	0.211417
2	3.90	0.000000
3	16.90	0.000000
4	73.15	0.000000
5	208.33	0.000000

Finally, to get the relative growth rate, we divide by $\lambda_0$.

lam_MG['relative growth rate'] = lam_MG['growth rate (1/hr)'] / lam0['MG']

# Take a look
lam_MG

	conc	growth rate (1/hr)	relative growth rate
0	0.21	0.405868	0.821553
1	0.90	0.211417	0.427947
2	3.90	0.000000	0.000000
3	16.90	0.000000	0.000000
4	73.15	0.000000	0.000000
5	208.33	0.000000	0.000000

We can do similar analysis for all of the other strains to get our growth rates.

Theoretical model for growth rate as a function of abx concentration

In the handout associated with the experiment, we worked out a theoretical model for growth rate as a function of external antibiotic concentration $a_\mathrm{ex}$. To solve for the growth rate, we have to find the solution to the dimensionless equation

$$\begin{align} 0 = \;&(K + \tilde{q}_\mathrm{ex, 0}^\mathrm{max} - 1)\tilde{\lambda}^3 \nonumber \\ &+ (2 - (1-K)(\tilde{P}_\mathrm{in} + \tilde{P}_\mathrm{out}) - (1-\kappa)(\tilde{q}_\mathrm{ex, 0}^\mathrm{max} + K)) \tilde{\lambda}^2 \nonumber \\ &-(1+\kappa(\tilde{q}_\mathrm{ex, 0}^\mathrm{max} + K) - (2 - (1-\kappa)K)\tilde{P}_\mathrm{out}-(1 - (1-2K)\kappa)\tilde{P}_\mathrm{in}\,\tilde{a}_\mathrm{ex})\tilde{\lambda} \nonumber \\ &+ (1+\kappa K)(\kappa \,\tilde{P}_\mathrm{in} \tilde{a}_\mathrm{ex} - \tilde{P}_\mathrm{out}), \end{align}$$

as derived in the handout. Upon finding the solution, we can get the dimensional growth rate by multiplying by $\lambda_0$.

By playing with parameters, the above cubic equation may have zero, one, two, or three roots on the interval $0 \le \lambda \le 1$. Though I will not derive it here, we can show that the following are true.

• Zero roots: In this situation there is no growth and bacteria die.
• One root: The value of $\lambda$ corresponding to this root represents a stable growth rate.
• Two roots: The greater of the two roots is stable and the lesser is unstable. There can be a population of cells growing at a rate given by the greater of the two roots and another population that is not growing.
• Three roots: The minimal and maximal roots are stable; the middle one is unstable. There is a population of cells growing at the faster rate and one growing at a slower rate.

In the case of two or three roots, the population growing at the faster rate will overwhelm the slow growers and will be the population we can see growing in the experiment. Therefore, our task is to find the largest root of the above cubic equation that is less than one. If no such root exists, the growth rate is zero.

Fortunately, Numpy has a built in function for finding roots of polynomials, np.polynomial.Polynomial.roots(). We can use this to find all three roots of the above cubic equation. The functions below accomplish this. Note that the default values for parameters known from the literature are included in the growth_rate() function.

def solve_growth_rate_dimensionless(
    aex, kappa, K, qmax, Pin, Pout, return_all_roots=False
):
    """Solve for the growth rate for a single antibiotic
    concentation in dimensionless units. All parameters
    are dimensionless. If `return_all_roots` is True,
    returns all physically viable steady state growth rates.
    Otherwise, only the single stable growth rate is returned.
    """
    # No antibiotic: max growth rate
    if aex == 0.0:
        lam = np.array([1.0])
    else:
        # Coefficients for cubic polynomial
        a3 = K + qmax - 1
        a2 = 2 - (1 - K) * (Pin * aex + Pout) - (1 - kappa) * (qmax + K)
        a1 = (
            -1
            - kappa * (qmax + K)
            + (2 - (1 - kappa) * K) * Pout
            + (1 - (1 - 2 * K) * kappa) * Pin * aex
        )
        a0 = (1 + kappa * K) * (kappa * Pin * aex - Pout)

        # Set up the cubic polynomial
        cubic = np.polynomial.Polynomial([a0, a1, a2, a3])

        # Solve the cubic
        lam = np.sort([r for r in cubic.roots() if np.isreal(r) and 0 <= r <= 1])

    # No roots, return 0
    if len(lam) == 0:
        return 0.0

    if return_all_roots:
        return lam if len(lam) > 1 else lam[0]
    else:
        return lam[-1]


def growth_rate(
    aex,            # µM
    qmax,           # µM/hr
    lam0,           # 1/hr
    rmin=19.3,      # µM
    kappa_t=0.061,  # 1/µM
    Pin=2.85,       # 1/hr
    Pout=2.85,      # 1/hr
    Kd=1.0,         # µM
    KM=10.0,        # µM
):
    """Solve for the growth rate for antibiotic concentration, or
    array of antibiotic concentrations, aex. Assumes all
    concentration units are µM and all time units are hours.
    """
    # Nondimensionalize
    kappa = kappa_t * rmin / lam0
    K = KM / Kd
    qmax = qmax / Kd / lam0
    Pin = Pin / lam0
    Pout = Pout / lam0
    aex = aex / Kd

    if np.isscalar(aex):
        result = solve_growth_rate_dimensionless(
            aex, kappa, K, qmax, Pin, Pout, return_all_roots=False
        )
    else:
        result = [
            solve_growth_rate_dimensionless(
                aex_val, kappa, K, qmax, Pin, Pout, return_all_roots=False
            )
            for aex_val in aex
        ]
        result = np.array(result)

    # Return dimensional growth rate
    return result * lam0

As an example, we can compute the growth rate for $q_\mathrm{ex, 0}^\mathrm{max} = 500$ with $\lambda_0 = 0.5$.

aex = np.linspace(0, 2, 2000)
lam = growth_rate(aex, 500, 0.5)

p = bokeh.plotting.figure(
    frame_width=350,
    frame_height=200,
    x_axis_label='aex (µM)',
    y_axis_label='λ  (1/hr)',
    x_range=[0, 2],
)

p.line(aex, lam, line_width=2)

bokeh.io.show(p)

The growth_rate() function is useful when doing a regression to find estimates for $q_{\mathrm{ex},0}^\mathrm{max}$ for each stain using the $\lambda/\lambda_0$ versus $a_\mathrm{ex}$ curves you obtain in this experiment. The function you will use in scipy.optimize.curve_fit() is shown below.

def fit_function(aex, qmax):
    """lam0 has to be externally defined as the abx-free growth rate."""
    return growth_rate(aex, qmax, lam0)