Drug combination optimization¶
NCI-ALMANAC · Cancer cell assays
In this tutorial, we’ll use GraphFLA to explore how drug combinations and doses affect cancer-cell growth. To bring categorical and ordered choices together, we’ll build a landscape from NCI-ALMANAC and examine the responses to different partners and dose levels.
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 ["almanac.csv", "almanac_doses.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
from IPython.display import display
2. Loading the dataset¶
NCI-ALMANAC is a large experimental screen of anticancer drug pairs. Researchers exposed cancer cell lines to combinations at different doses and measured cell growth. To explore one manageable search space, we’ll focus on 786-0 cells treated with idarubicin hydrochloride and one of 100 partner drugs, using measurements from screener 1A.
Each drug has three tested dose levels, giving 100 × 3 × 3 = 900 configurations. We’ll minimise time-zero-corrected PercentGrowth, so smaller values mean less growth in the assay. This is a growth response, not a synergy score.
| Column | Meaning |
|---|---|
partner_nsc |
Partner-drug identifier. |
anchor_dose_index, partner_dose_index |
Ordered low, middle and high tested doses. |
mean_percent_growth |
Original PercentGrowth response. |
Let’s load this complete block. Each point has one observation despite the mean_ prefix. The CSV retains identifiers for the 12 studies and 300 plates it spans.
df = pd.read_csv("data/almanac.csv")
df[["partner_name", "anchor_dose_index", "partner_dose_index", "mean_percent_growth"]].head()
Output
| partner_name | anchor_dose_index | partner_dose_index | mean_percent_growth | |
|---|---|---|---|---|
| 0 | Methotrexate | 1 | 1 | 67.971 |
| 1 | Methotrexate | 1 | 2 | 17.190 |
| 2 | Methotrexate | 1 | 3 | 16.230 |
| 3 | Methotrexate | 2 | 1 | 81.874 |
| 4 | Methotrexate | 2 | 2 | 15.997 |
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 partner identity and two dose indices form X; PercentGrowth forms f. Anchor doses are 0.1, 1 and 10 nM. Partner concentrations depend on drug identity, so an index does not represent a universal concentration. The accompanying lookup preserves the actual doses.
dose_lookup = pd.read_csv("data/almanac_doses.csv")
display(dose_lookup.head())
X = df[["partner_nsc", "anchor_dose_index", "partner_dose_index"]]
f = df["mean_percent_growth"]
X.head()
Output
| component_role | compound_nsc | compound_name | dose_index | concentration | unit | |
|---|---|---|---|---|---|---|
| 0 | anchor | 256439 | Idarubicin hydrochloride | 1 | 1.000000e-10 | M |
| 1 | anchor | 256439 | Idarubicin hydrochloride | 2 | 1.000000e-09 | M |
| 2 | anchor | 256439 | Idarubicin hydrochloride | 3 | 1.000000e-08 | M |
| 3 | partner | 740 | Methotrexate | 1 | 1.000000e-08 | M |
| 4 | partner | 740 | Methotrexate | 2 | 1.000000e-07 | M |
| partner_nsc | anchor_dose_index | partner_dose_index | |
|---|---|---|---|
| 0 | 740 | 1 | 1 |
| 1 | 740 | 1 | 2 |
| 2 | 740 | 1 | 3 |
| 3 | 740 | 2 | 1 |
| 4 | 740 | 2 | 2 |
4. Constructing the landscape¶
We can use Landscape with one categorical variable and two ordinal variables. Neighbours replace the partner at fixed dose indices or move one dose to its next sampled level. Switching partners may also change its physical concentration.
We set maximize=False because smaller growth responses are better for this objective. Cell line and anchor identity remain fixed.
landscape = Landscape(maximize=False)
landscape.build_from_data(
X, f, data_types={
"partner_nsc": "categorical",
"anchor_dose_index": "ordinal",
"partner_dose_index": "ordinal",
}, verbose=False,
)
Output
Landscape(kind='default', maximize=False)
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, 900 configurations, 45746 edges, 3 local optima
For a closer look at individual drug-and-dose 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
| partner_nsc | anchor_dose_index | partner_dose_index | fitness | plateau_id | plateau_size | in_degree | out_degree | is_lo | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 740 | 1 | 1 | 67.971 | -1 | 1 | 98 | 3 | False |
| 1 | 740 | 1 | 2 | 17.190 | -1 | 1 | 99 | 3 | False |
| 2 | 740 | 1 | 3 | 16.230 | -1 | 1 | 85 | 16 | False |
| 3 | 740 | 2 | 1 | 81.874 | -1 | 1 | 98 | 4 | False |
| 4 | 740 | 2 | 2 | 15.997 | -1 | 1 | 102 | 1 | False |
We can also summarise the growth responses with pandas:
landscape.get_data()["fitness"].describe()
Output
count 900.000000
mean 78.274919
std 30.652225
min -65.403000
25% 66.784250
50% 89.030500
75% 100.652500
max 118.257000
Name: fitness, dtype: float64
To inspect the combination with the lowest recorded growth response, we can use the global-optimum node index:
landscape[landscape.go_index]
Output
{'partner_nsc': np.int64(34462),
'anchor_dose_index': np.int64(3),
'partner_dose_index': np.int64(3),
'fitness': np.float64(-65.403),
'plateau_id': -1,
'plateau_size': 1,
'in_degree': 101,
'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.2737
6.3 Walsh–Hadamard decomposition¶
Does changing the dose have the same effect with every partner? walsh_hadamard() fits individual-variable effects and pairwise interactions, treating each tested dose level as a discrete state. delta_r2 shows how much training variation each order adds.
Three-way effects remain outside this fit. These interactions describe the selected response table; they are not pharmacological synergy estimates. The model-variance spectrum describes uniform combinations within the fitted model, separately from cumulative training R².
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 | 30.635191 | 1 | 1 | NaN | 0.000000 | <NA> |
| 1 | 1 | 0.658349 | 0.658349 | 17.906563 | 104 | 104 | NaN | 0.664854 | <NA> |
| 2 | 2 | 0.990215 | 0.331866 | 3.030409 | 504 | 504 | NaN | 0.335146 | <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 | |
|---|---|---|---|---|
| 247 | 2 | (1, 3) | 740_1_34462-1_3_3 | -97.589333 |
| 403 | 2 | (1, 3) | 740_1_747599-1_3_3 | 72.011000 |
| 402 | 2 | (1, 3) | 740_1_747599-1_3_2 | 71.959000 |
| 127 | 2 | (1, 3) | 740_1_122758-1_3_3 | 66.875333 |
| 499 | 2 | (1, 3) | 740_1_92859-1_3_3 | 65.705667 |
| 334 | 2 | (1, 3) | 740_1_702294-1_3_2 | 65.179667 |
| 395 | 2 | (1, 3) | 740_1_743414-1_3_3 | 64.827000 |
| 150 | 2 | (1, 3) | 740_1_13875-1_3_2 | 63.879333 |
6.4 Fitness-distance correlation¶
Distance adds partner replacement to the number of sampled dose-rank steps. It does not compare absolute molar concentrations across different drugs.
We use Spearman fdc to compare ranks. A positive value means lower 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.5809
Pairwise terms increase training R² from 0.658 to 0.990. This captures much of the joint partner-and-dose response variation, while the positive FDC (0.581) associates proximity to the selected optimum with lower growth.
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.273671
fdc 0.580901
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¶
Holbeck S. L. et al. (2017). The National Cancer Institute ALMANAC: A Comprehensive Screening Resource for the Detection of Anticancer Drug Pairs with Enhanced Therapeutic Activity. We use a complete fixed-context block from the official NCI-ALMANAC data, with source study, plate and concentration fields retained.