Chemical reaction optimization¶
Suzuki–Miyaura coupling
In this tutorial, we’ll use GraphFLA to explore how reaction conditions affect the outcome of a chemical reaction. To see how this works, we’ll build and analyse a landscape from a Suzuki–Miyaura reaction screen.
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 ["suzuki.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¶
Suzuki–Miyaura coupling joins two molecular fragments by forming a carbon–carbon bond. Finding suitable reaction conditions can involve trying many combinations of reagents and solvents. To explore this search problem, we’ll use the dataset from Perera et al. (2018), who developed an automated flow platform and screened 5,760 reactions using small amounts of material.
We’ll focus on one pair of starting materials, 6-chloroquinoline and boronic acid 2a, so we can compare conditions for making the same product. This gives 12 ligand choices × 8 base/reagent choices × 4 solvent labels = 384 reactions. Our objective will be to maximise the reported UV product-area percentage, a screening readout rather than an isolated product yield.
| Column | Meaning |
|---|---|
ligand |
Ligand used with the catalyst, including a no-addition control. |
base |
Base/reagent identity, including a no-addition control. |
solvent |
Solvent condition. |
response_uv_pct |
Product-area percentage measured by UV. |
Let’s load these measurements. None denotes a no-addition control, so we’ll preserve it as a category.
df = pd.read_csv("data/suzuki.csv", keep_default_na=False, float_precision="round_trip")
df.head()
Output
| Reaction_No | reactant_1 | reactant_2 | ligand | base | solvent | response_uv_pct | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | 6-chloroquinoline | 2a, Boronic Acid | P(tBu)3 | NaOH | MeCN | 4.764109 |
| 1 | 2 | 6-chloroquinoline | 2a, Boronic Acid | P(Ph)3 | NaOH | MeCN | 4.120962 |
| 2 | 3 | 6-chloroquinoline | 2a, Boronic Acid | AmPhos | NaOH | MeCN | 2.583837 |
| 3 | 4 | 6-chloroquinoline | 2a, Boronic Acid | P(Cy)3 | NaOH | MeCN | 4.443171 |
| 4 | 5 | 6-chloroquinoline | 2a, Boronic Acid | P(o-Tol)3 | NaOH | MeCN | 1.949874 |
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.
The three condition identities form X, and the UV response forms f. The substrates remain fixed, so they are not additional variables in this landscape.
X = df[["ligand", "base", "solvent"]]
f = df["response_uv_pct"]
X.head()
Output
| ligand | base | solvent | |
|---|---|---|---|
| 0 | P(tBu)3 | NaOH | MeCN |
| 1 | P(Ph)3 | NaOH | MeCN |
| 2 | AmPhos | NaOH | MeCN |
| 3 | P(Cy)3 | NaOH | MeCN |
| 4 | P(o-Tol)3 | NaOH | MeCN |
4. Constructing the landscape¶
We can use Landscape with three categorical variables. Neighbouring reactions replace one ligand, base or solvent while retaining the other two choices. A category change does not imply chemical similarity between the alternatives.
We set maximize=True because a larger UV product-area response 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'): 3 variables, 384 configurations, 4020 edges, 4 local optima
For a closer look at individual reaction conditions, 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
| ligand | base | solvent | fitness | plateau_id | plateau_size | in_degree | out_degree | is_lo | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | P(tBu)3 | NaOH | MeCN | 4.764109 | -1 | 1 | 9 | 12 | False |
| 1 | P(Ph)3 | NaOH | MeCN | 4.120962 | -1 | 1 | 4 | 17 | False |
| 2 | AmPhos | NaOH | MeCN | 2.583837 | -1 | 1 | 8 | 13 | False |
| 3 | P(Cy)3 | NaOH | MeCN | 4.443171 | -1 | 1 | 5 | 16 | False |
| 4 | P(o-Tol)3 | NaOH | MeCN | 1.949874 | -1 | 1 | 1 | 20 | False |
We can also summarise the UV product-area responses with pandas:
landscape.get_data()["fitness"].describe()
Output
count 384.000000
mean 23.922465
std 24.905576
min 0.000000
25% 5.835257
50% 14.149675
75% 32.757138
max 93.674515
Name: fitness, dtype: float64
To inspect the reaction with the highest measured UV response, we can use the global-optimum node index:
landscape[landscape.go_index]
Output
{'ligand': 'XPhos',
'base': 'NaOH',
'solvent': 'MeCN',
'fitness': np.float64(93.6745145376158),
'plateau_id': -1,
'plateau_size': 1,
'in_degree': 21,
'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: 4
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.3979
6.3 Walsh–Hadamard decomposition¶
Do the three choices contribute independently, or does their combination matter? walsh_hadamard() fits additive terms and interactions between pairs of variables. Compare the first- and second-order rows: delta_r2 shows the extra training variation explained by pairwise terms.
The fit stops at order two, leaving three-way interactions unresolved. model_variance_fraction describes the fitted model’s variance under uniform combinations; it is a different quantity from delta_r2.
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.000000 | 0.000000 | 24.873126 | 1 | 1 | NaN | 0.000000 | <NA> |
| 1 | 1 | 0.636035 | 0.636035 | 15.005832 | 22 | 22 | NaN | 0.662877 | <NA> |
| 2 | 2 | 0.959507 | 0.323472 | 5.005190 | 153 | 153 | NaN | 0.337123 | <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 | |
|---|---|---|---|---|
| 104 | 2 | (1, 3) | P(tBu)3_1_SPhos-MeCN_3_MeOH | -67.331393 |
| 114 | 2 | (1, 3) | P(tBu)3_1_XPhos-MeCN_3_MeOH | -65.901641 |
| 115 | 2 | (1, 3) | P(tBu)3_1_XPhos-MeCN_3_THF | -55.006236 |
| 65 | 2 | (1, 3) | P(tBu)3_1_None-MeCN_3_THF | -47.036022 |
| 125 | 2 | (1, 3) | P(tBu)3_1_Xantphos-MeCN_3_THF | -46.835056 |
| 95 | 2 | (1, 3) | P(tBu)3_1_P(o%2DTol)3-MeCN_3_THF | -46.590488 |
| 105 | 2 | (1, 3) | P(tBu)3_1_SPhos-MeCN_3_THF | -45.535031 |
| 85 | 2 | (1, 3) | P(tBu)3_1_P(Ph)3-MeCN_3_THF | -45.014717 |
6.4 Fitness-distance correlation¶
Distance counts how many condition choices differ from the best observed reaction. It does not measure molecular 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.0535
Here, the pairwise fit reaches R² ≈ 0.960, but FDC is close to zero (0.054). Explaining response variation does not necessarily give a simple distance-based guide towards the best reaction.
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.397925
fdc 0.053540
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¶
Perera D. et al. (2018). A platform for automated nanomole-scale reaction screening and micromole-scale synthesis in flow. We use one complete substrate-pair grid from Data File S1.