Enzyme inhibitor discovery¶

Cyanimide library · Mouse USP18

In this tutorial, we’ll use GraphFLA to explore how combining chemical building blocks affects enzyme inhibition. To follow these choices from molecules to activity, we’ll build and analyse a landscape from a library of cyanimide-containing compounds.

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 ["cyanimide.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¶

One way to search for enzyme inhibitors is to assemble a library of compounds from interchangeable building blocks, then measure their activity. To explore this approach, we’ll use the screen from Kooij et al. (2025). The researchers coupled 16 cyanimide-bearing amines with 469 carboxylic acids and tested the resulting reaction mixtures against several proteases, enzymes that cleave proteins.

We’ll focus on mouse USP18 and maximise its inhibition percentage at the screening concentration of 1.25 µM. The measurements describe crude reaction mixtures, rather than purified compounds’ dose–response potency.

Column	Meaning
`bb`, `ca`	Amine and carboxylic-acid building-block IDs.
`mUSP18_mean_percent_inhibition`	Mean primary-screen inhibition.
`n_wells`	Number of assay wells for that combination.

Let’s load the 7,504 building-block combinations. The table averages the 32 pairs with duplicate wells and records the number of wells for each pair.

df = pd.read_csv("data/cyanimide.csv")
df[["bb", "ca", "mUSP18_mean_percent_inhibition", "n_wells"]].head()

Output

	bb	ca	mUSP18_mean_percent_inhibition	n_wells
0	BB01	CA001	21.663005	1
1	BB01	CA002	21.622812	1
2	BB01	CA003	22.786787	1
3	BB01	CA004	26.376782	1
4	BB01	CA005	20.430703	1

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 of X.

The building-block identities form X, and mean mUSP18 inhibition forms f. Each pair now occurs once. We retain negative readings near the assay baseline and use the original inhibition scale.

X = df[["bb", "ca"]]
f = df["mUSP18_mean_percent_inhibition"]
X.head()

Output

	bb	ca
0	BB01	CA001
1	BB01	CA002
2	BB01	CA003
3	BB01	CA004
4	BB01	CA005

4. Constructing the landscape¶

We can use Landscape with two categorical variables. Neighbouring compounds share one building block and replace the other. The IDs identify reagents; their ordering does not encode chemical distance.

We set maximize=True to favour stronger inhibition. Selectivity across the protease panel remains a separate 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, 7504 configurations, 1812150 edges, 6 local optima

For a closer look at individual building-block 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

	bb	ca	fitness	plateau_id	plateau_size	in_degree	out_degree	is_lo
0	BB01	CA001	21.663005	-1	1	226	257	False
1	BB01	CA002	21.622812	-1	1	227	256	False
2	BB01	CA003	22.786787	-1	1	235	248	False
3	BB01	CA004	26.376782	-1	1	264	219	False
4	BB01	CA005	20.430703	-1	1	221	262	False

We can also summarise the inhibition measurements with pandas:

landscape.get_data()["fitness"].describe()

Output

count    7504.000000
mean       15.545357
std        19.983758
min        -5.618241
25%         2.926512
50%         8.092254
75%        20.256774
max        99.172273
Name: fitness, dtype: float64

To inspect the combination with the highest mean inhibition, we can use the global-optimum node index:

landscape[landscape.go_index]

Output

{'bb': 'BB01',
 'ca': 'CA148',
 'fitness': np.float64(99.17227342848624),
 'plateau_id': -1,
 'plateau_size': 1,
 'in_degree': 483,
 '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: 6

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.3418

6.3 Walsh–Hadamard decomposition¶

Can independent building-block contributions explain the screen? We first use walsh_hadamard() at order one across all 7,504 compounds. Its training r2 measures the additive fit; the remaining variation is not assigned to specific interaction coefficients.

A full second-order fit needs 7,504 coefficients. To demonstrate pairwise decomposition at a manageable size, we also use all 16 amines with the first 12 acid IDs, selected without looking at inhibition. Its order summary describes only this 192-compound slice.

wh = analysis.walsh_hadamard(landscape, max_order=1)
wh["order_summary"]

Output

	order	r2	delta_r2	rmse	n_terms	rank	alpha	model_variance_fraction	n_nonzero
0	0	0.000000	0.000000	19.982426	1	1	NaN	0.0	<NA>
1	1	0.679878	0.679878	11.305923	484	484	NaN	1.0	<NA>

acid_ids = sorted(df["ca"].unique())[:12]
slice_df = df[df["ca"].isin(acid_ids)]
wh_landscape = Landscape(maximize=True).build_from_data(
    slice_df[["bb", "ca"]], slice_df["mUSP18_mean_percent_inhibition"],
    data_types={"bb": "categorical", "ca": "categorical"}, verbose=False,
)
wh_slice = analysis.walsh_hadamard(wh_landscape, max_order=2)
print(f"Pairwise example: {wh_landscape.n_configs} compounds; acid IDs {acid_ids}")
wh_slice["order_summary"]

Output

Pairwise example: 192 compounds; acid IDs ['CA001', 'CA002', 'CA003', 'CA004', 'CA005', 'CA006', 'CA007', 'CA008', 'CA009', 'CA010', 'CA011', 'CA012']

	order	r2	delta_r2	rmse	n_terms	rank	alpha	model_variance_fraction	n_nonzero
0	0	0.000000	0.000000	1.386260e+01	1	1	NaN	0.000000	<NA>
1	1	0.802455	0.802455	6.161371e+00	27	27	NaN	0.802455	<NA>
2	2	1.000000	0.197545	4.012568e-14	192	192	NaN	0.197545	<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
49	1	(2,)	CA001_2_CA035	58.010441
82	1	(2,)	CA001_2_CA068	52.601033
421	1	(2,)	CA001_2_CA407	44.520320
206	1	(2,)	CA001_2_CA192	43.728291
75	1	(2,)	CA001_2_CA061	42.928205
263	1	(2,)	CA001_2_CA249	42.319466
79	1	(2,)	CA001_2_CA065	40.645378
306	1	(2,)	CA001_2_CA292	39.326892

6.4 Fitness-distance correlation¶

Distance counts replaced building blocks relative to the best observed pair. It does not measure molecular fingerprint similarity. Local optima, autocorrelation and FDC all use the full library.

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.2367

The full-library additive fit explains about 68.0% of variation, with six local optima. The 19.8% pairwise contribution in the smaller example applies only to that selected slice; it is not the whole-library interaction fraction.

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.341791
fdc               -0.236737
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¶

Kooij R. et al. (2025). High-Throughput Synthesis and Screening of a Cyanimide Library Identifies Selective Inhibitors of ISG15-Specific Protease mUSP18. We use the supplementary primary-screen measurements, averaged by building-block pair.