Microbial growth optimization¶
BacPUS · Gut bacteria and polysaccharides
In this tutorial, we’ll use GraphFLA to explore how nutrient choices affect bacterial growth. To compare these responses, we’ll build a landscape from gut bacterial isolates grown on different polysaccharides, then examine which strain–substrate combinations support stronger growth.
Download notebook Notebook + data
Run the notebook in Google Colab, or download it to run locally. It downloads its data when no data/ folder is beside it.
1. Setting up¶
On Google Colab, the cell below installs GraphFLA. In any environment, it also downloads the dataset into a data/ folder next to the notebook if the files are not already there.
We then import GraphFLA for landscape construction and analysis, and pandas for working with the data.
import sys
from pathlib import Path
from urllib.request import urlretrieve
if "google.colab" in sys.modules:
%pip install -q graphfla==0.4.0
DATA_URL = "https://raw.githubusercontent.com/COLA-Laboratory/GraphFLA/v0.4.0/tutorials/datasets/data/"
Path("data").mkdir(exist_ok=True)
for name in ["bacpus.csv"]:
if not Path("data", name).exists():
urlretrieve(DATA_URL + name, Path("data", name))
from graphfla import analysis
from graphfla.landscape import Landscape
import pandas as pd
2. Loading the dataset¶
Gut bacteria differ in their ability to grow on polysaccharides, carbohydrates made from chains of sugar units. To examine how growth depends on both the bacterium and its nutrient source, we’ll use the BacPUS screen from Qu et al. (2025). The researchers cultured 28 isolates separately on 20 polysaccharide preparations, measuring all 560 combinations.
We’ll use mean optical density at 600 nm after 48 hours (OD600) as the growth readout, with larger values as our objective. The table averages available replicates and retains small negative readings after background correction.
| Column | Meaning |
|---|---|
strain_id |
Identifier distinguishing individual bacterial isolates. |
polysaccharide |
Identifier of the tested preparation. |
OD600_48h |
Mean growth readout after 48 hours. |
Let’s load the measurements, keeping in mind that each response comes from one isolate grown on one substrate.
df = pd.read_csv("data/bacpus.csv")
df[["strain_id", "strain", "polysaccharide", "OD600_48h"]].head()
Output
| strain_id | strain | polysaccharide | OD600_48h | |
|---|---|---|---|---|
| 0 | DA103 | Bacteroides_massiliensis | APs | 0.030000 |
| 1 | DA103 | Bacteroides_massiliensis | AdPs | -0.001833 |
| 2 | DA103 | Bacteroides_massiliensis | AnPs | 0.010833 |
| 3 | DA103 | Bacteroides_massiliensis | CPPs | -0.003333 |
| 4 | DA103 | Bacteroides_massiliensis | DPs_1 | 0.084500 |
3. Preparing the inputs¶
To construct a landscape, GraphFLA needs two aligned inputs:
X: one row per configuration and one column per variable.f: one measured or calculated outcome for each row ofX.
We select the isolate and polysaccharide identifiers as X, and OD600 as f. Strain IDs keep distinct isolates of the same species separate; species names alone would merge different assay configurations.
X = df[["strain_id", "polysaccharide"]]
f = df["OD600_48h"]
X.head()
Output
| strain_id | polysaccharide | |
|---|---|---|
| 0 | DA103 | APs |
| 1 | DA103 | AdPs |
| 2 | DA103 | AnPs |
| 3 | DA103 | CPPs |
| 4 | DA103 | DPs_1 |
4. Constructing the landscape¶
We can use Landscape with two categorical variables. Neighbours replace the strain or the polysaccharide while keeping the other choice fixed. These moves compare assay configurations; they do not represent mutations or evolutionary transitions between species.
We set maximize=True because greater growth is our objective.
landscape = Landscape(maximize=True)
landscape.build_from_data(
X, f, data_types={column: "categorical" for column in X},
neighborhood_strategy="active", verbose=False,
)
Output
Landscape(kind='default', maximize=True)
5. Inspecting the landscape¶
Let’s first look at what we’ve built. Printing the landscape shows its variables, configurations, improving edges and local optima:
print(landscape)
Output
Landscape(kind='default'): 2 variables, 560 configurations, 12832 edges, 3 local optima
For a closer look at individual strain–substrate combinations, we can call get_data(). The outcome appears as fitness; out_degree counts directly improving moves, and is_lo indicates local-optimum membership.
landscape.get_data().head()
Output
| strain_id | polysaccharide | fitness | plateau_id | plateau_size | in_degree | out_degree | is_lo | |
|---|---|---|---|---|---|---|---|---|
| 0 | DA103 | APs | 0.030000 | -1 | 1 | 14 | 32 | False |
| 1 | DA103 | AdPs | -0.001833 | -1 | 1 | 6 | 40 | False |
| 2 | DA103 | AnPs | 0.010833 | -1 | 1 | 11 | 35 | False |
| 3 | DA103 | CPPs | -0.003333 | -1 | 1 | 3 | 43 | False |
| 4 | DA103 | DPs_1 | 0.084500 | -1 | 1 | 35 | 11 | False |
We can also summarise the growth measurements with pandas:
landscape.get_data()["fitness"].describe()
Output
count 560.000000
mean 0.132891
std 0.146531
min -0.011667
25% 0.020167
50% 0.074667
75% 0.191250
max 0.692667
Name: fitness, dtype: float64
To inspect the strain–substrate combination with the highest growth score, we can use the global-optimum node index:
landscape[landscape.go_index]
Output
{'strain_id': 'DA320',
'polysaccharide': 'GPs',
'fitness': np.float64(0.6926666666666668),
'plateau_id': -1,
'plateau_size': 1,
'in_degree': 46,
'out_degree': 0,
'is_lo': True}
6. Analysing the landscape¶
Next, we’ll count local optima, examine local fitness similarity and interaction orders, and assess the trend towards the best observed configuration.
6.1 Number of local optima¶
We can start with landscape.n_lo. Each connected neutral optimum plateau counts once. A local optimum has no improving move out of it, including through its neutral plateau. Multiple optima mean successive improvements can end at different configurations.
print(f"Number of local optima: {landscape.n_lo}")
Output
Number of local optima: 3
6.2 Fitness autocorrelation¶
How similar are responses at neighbouring steps? autocorrelation() measures persistence along sampled walks. Higher values indicate more persistence under the same walk settings.
We use 200 walks of up to 20 visited configurations, lag one and a fixed seed. Walks traverse stored edges in either direction; separately stored neutral pairs are excluded.
walk_autocorrelation = analysis.autocorrelation(
landscape, walk_length=20, walk_times=200, lag=1, seed=42,
)
print(f"Lag-1 fitness autocorrelation: {walk_autocorrelation:.4f}")
Output
Lag-1 fitness autocorrelation: 0.1984
6.3 Walsh–Hadamard decomposition¶
A polysaccharide may support one isolate much better than another. walsh_hadamard() separates independent strain and substrate contributions from their pairwise interaction. Compare the first-order fit with the second-order fit using r2 and delta_r2.
With two variables and a complete table, the second-order model can reproduce every measured value. This exact fit is descriptive, not evidence of predictive accuracy or interactions between organisms in a community. The uniform-product spectrum summarises variance within the fitted model.
wh = analysis.walsh_hadamard(landscape, max_order=2)
wh["order_summary"]
Output
| order | r2 | delta_r2 | rmse | n_terms | rank | alpha | model_variance_fraction | n_nonzero | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0.00000 | 0.00000 | 1.464001e-01 | 1 | 1 | NaN | 0.00000 | <NA> |
| 1 | 1 | 0.50217 | 0.50217 | 1.032956e-01 | 47 | 47 | NaN | 0.50217 | <NA> |
| 2 | 2 | 1.00000 | 0.49783 | 1.251036e-15 | 560 | 560 | NaN | 0.49783 | <NA> |
We can also inspect the largest fitted nonconstant coefficients. Their units follow the response, and their signs depend on the encoded contrasts:
wh["coefficients"].query("order > 0").sort_values(
"coefficient", key=abs, ascending=False,
).head(8)
Output
| order | positions | term | coefficient | |
|---|---|---|---|---|
| 320 | 2 | (1, 2) | DA103_1_DA320-APs_2_GPs | 0.581833 |
| 369 | 2 | (1, 2) | DA103_1_DA486-APs_2_SPs | 0.491000 |
| 104 | 2 | (1, 2) | DA103_1_DA1247-APs_2_AdPs | 0.480167 |
| 217 | 2 | (1, 2) | DA103_1_DA1479-APs_2_SPs | 0.464500 |
| 403 | 2 | (1, 2) | DA103_1_DA557-APs_2_PSPs | 0.443500 |
| 408 | 2 | (1, 2) | DA103_1_DA57-APs_2_AdPs | 0.432000 |
| 434 | 2 | (1, 2) | DA103_1_DA647-APs_2_GPs | 0.422333 |
| 389 | 2 | (1, 2) | DA103_1_DA557-APs_2_AdPs | 0.419167 |
6.4 Fitness-distance correlation¶
Distance counts whether strain and substrate choices differ from the best observed combination. It does not encode taxonomic or chemical similarity.
We use Spearman fdc to compare ranks. A negative value means higher responses tend to occur nearer the optimum. A value near zero indicates little monotonic association. This trend does not guarantee an improving path from every configuration.
fitness_distance_r = analysis.fdc(landscape, method="spearman")
print(f"Fitness-distance correlation: {fitness_distance_r:.4f}")
Output
Fitness-distance correlation: -0.2184
Independent strain and substrate terms explain about 50.2% of variation. Pairwise terms describe the remaining 49.8% in this complete table, showing substantial dependence of growth on the particular strain–substrate combination.
7. Running several analyses together¶
We can collect autocorrelation and FDC with analysis.profile(), keeping the same walk settings. landscape.n_lo supplies the count directly; the W–H result keeps its separate coefficient and order-summary tables.
analysis.profile(
landscape,
metrics=["autocorrelation", "fdc"],
params={"autocorrelation": {"walk_length": 20, "walk_times": 200, "lag": 1}},
seed=42, progress=False,
)
Output
autocorrelation 0.198431
fdc -0.218363
dtype: float64
8. Analysis reference¶
For further exploration, analysis.list_metrics() lists the metrics available to profile(). The worked W–H result also exposes coefficients and fit_info; its order summary separates cumulative training fit from the fitted model’s variance spectrum.
Data sources¶
Qu et al. (2025). BacPUS strain–polysaccharide study. The CSV contains isolate identifiers and mean 48-hour OD600 measurements.