Screening too many chemicals to count

Searching trillions of structures on a laptop

How to predict on chemical space without enumerating it

In modern AI-aided drug discovery, the chemical space addressed by the drug discovery funnel is truly vast. Virtual screening of large libraries of compounds is much higher-throughput than traditional primary screening. However, the chemical space that can be explored is still limited by the ability to enumerate and score compounds; this does not reasonably scale over many billions of compounds because of the expense of calculation on each compound, no matter how fast the calculations are.

To address this problem, we exploit the fact that we only want the highest scoring structures from a library. This is a very similar situation to online marketing or advertising matching, in which recommendations need to be made from relatively little data. There is a deep understanding of the methods that succeed in such a task. One prominent example is also one of the oldest and most versatile - Thompson sampling.

We’ll explore how Thompson sampling can be used to score chemical structures without enumerating them, and how it can be applied to the problem of virtual screening in drug discovery. For this application in chemistry, we’ll use a convenient implementation from a recent paper (Klarich et al. 2024), which is proving to be a valuable tool to exponentially amplify the power virtual of screening, by screening over fragments rather than over structures explicitly.

Defining a chemical space implicitly

We continue with the recurring example of ionizable lipids in lipid nanoparticles (LNPs) for drug delivery.

The celebrated Pfizer-BioNTech and Moderna mRNA vaccines for COVID-19 use modular lipid nanoparticles (LNPs) to deliver mRNA. These formulations use different ionizable lipids, known as ALC-0315 (Pfizer-BioNTech) and SM-102 (Moderna), which are often looked upon in a modular manner:

(a) Jorgensen et al. 2023

Figure 1: Adapted from Jorgensen et al. 2023.

Their “reaction skeletons” are similar (ester linkages connecting two tails with some branching to a tertiary amine head group) but have important differences (e.g. their ester groups face opposite directions with respect to the head group). Suppose we want to try varying these components from one of these basic structures.

Take SM-102 as an example (for technical reasons (Zhang et al. 2023)). We might want to vary the head and tails, keeping the ester linkers intact. Let’s say we want to modify it precisely, using the large and diverse collection of hundreds of primary amine heads and hydrophobic tails outlined in previous posts. We will perform the following modifications:

• Try a number of different head groups corresponding to different amines.

• Keep the alkyl chain spacers separating the acid groups from the amine nitrogen (currently each tail has 5 spacer carbons before the carbonyl group).

• Try different tail groups (the two alcohol-derived portions), varying the length, branching, and saturation.

LNPs are a versatile and modular technology in which the IL is formulated together with multiple (normally 3) other chemical components - they are engineered to work together to deliver the drug to its target. So the IL is not fully determinative. Such details are out of the scope of our presentation, and they may be more important for industrial-grade drug discovery efforts. However, the machine learning and tools presented here are generically useful in production-ready pipelines, and form the backbone of even more advanced engineering efforts.

In silico “reaction” rule

To implement Thompson sampling, we need to write these modifications as a multi-component in silico “reaction” in an unambiguous way. This will have several components per the above description:

• A head group, which is a primary amine.
• The first tail (as an alcohol).
• The second tail (as an alcohol).

We insert the alkyl spacer carbons on either side.

from rdkit import Chem
from rdkit.Chem import rdChemReactions


rxn_str = "[#6:1][NX3;H2:2].[#6:3][OX2H1:4].[#6:5][OX2H1:6]>>[#6:1][N:2](CCCCCCCC(=O)[OX2H0:4][#6:3])CCCCCC(=O)[OX2H0:6][#6:5]"
rdChemReactions.ReactionFromSmarts(rxn_str)

Tail groups

Next, we gather a set of tail groups we might want to try in this variant screening. Drug designers in this field vary the tail groups in a number of ways, and we include a representative sample of these variations:

• Varying the branching of the tail (unbranched (1-branched) and 2-branched tails explored)
• Varying the length of the tail (6-16 carbons as is standard, shorter for branched tails)
• Varying the saturation of the tail (each branch of each tail can contain a double bond somewhere along its length)

We can systematically generate a set of tail groups that vary in these ways. For virtual screening purposes, we’ll do this combinatorially, generating all possible combinations of these variations. This is discussed in more detail in another post.

In this case, the process results in thousands of candidate tail substructures.

source_pfx = "../../files/chem/"
sourcing_tools_path = source_pfx + "sourcing_ionizable_lipid_structures.py"

from importlib.machinery import SourceFileLoader
IL_structures = SourceFileLoader("IL_structures", sourcing_tools_path).load_module()

from itertools import combinations, product

from rdkit import RDLogger
lg = RDLogger.logger()
lg.setLevel(RDLogger.CRITICAL) # Set the logging level to CRITICAL to suppress all messages below this level


# design‑space parameters -----------------------------------------------------
unbranched_lengths        = [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]   # carbons / tail
branched_lengths        = [6, 7, 8, 9, 10, 11, 12]   # carbons / tail
unsat_choices  = (0, 1, 2)                 # number of C=C per tail
# -----------------------------------------------------------------------------

frags = set()


# unbranched tails
for L, d in product(unbranched_lengths, unsat_choices):
    interior = range(1, L - 1)                   # keep core‑C and C‑OH single
    if d > len(interior):
        continue

    for db in combinations(interior, d):
        if not IL_structures.spacing_ok(db):
            continue
        try:
            smi = IL_structures.build_linear_tail(L, db)
            Chem.SanitizeMol(Chem.MolFromSmiles(smi))
            frags.add(smi)
        except:
            pass                           # drops impossible ones silently

# ---------- branched ------------
for L1, L2, d1, d2 in product(branched_lengths, branched_lengths, unsat_choices, unsat_choices):
    interior1 = range(1, L1)                     # skip bond 0 (core‑C1)
    interior2 = range(1, L2)

    if d1 > len(interior1) or d2 > len(interior2):
        continue

    for db1 in combinations(interior1, d1):
        if not IL_structures.spacing_ok(db1):
            continue
        for db2 in combinations(interior2, d2):
            if not IL_structures.spacing_ok(db2):
                continue
            try:
                smi = IL_structures.build_branched_tail(L1, L2, db1, db2)
                Chem.SanitizeMol(Chem.MolFromSmiles(smi))
                frags.add(smi)
            except Exception:
                pass

We can write this to a file once we’re done enumerating these moieties.

import pandas as pd

tail_pfx = "../../files/LNP_catalog/"
tail_file_path = tail_pfx + "tail_frags.smi"
names_tails = []

frags_tails = [Chem.CanonSmiles(f) for f in frags if Chem.MolFromSmiles(f) is not None]
names_tails.extend(["tail-" + str(x) for x in range(len(frags_tails))])
tails_df = pd.DataFrame([frags_tails, names_tails], index=["SMILES", "NAME"]).T
tails_df.to_csv(tail_file_path, index=False, header=False, sep=" ")

print(f"{len(frags_tails)} unique fragments written.")

9892 unique fragments written.

Note that this type of combinatorial expansion is much easier to do computationally than in an assay, playing to the strengths of virtual screening.

The implementation involves helper functions linear_smiles and branched_smiles, which automatically generate tails of a particular nature. When all the combinations are expanded, we have a set of many hundreds of tail groups that can be used in the screening.

import mols2grid
mols2grid.display([Chem.MolFromSmiles(x) for x in frags_tails],mol_col="mol", n_cols=7, n_rows=4)

Head groups

The reaction skeleton we’ve defined above involves a primary amine reactant that becomes the head group of the IL. Due to the importance of the head group to the therapeutic behavior of the IL, it is often the main focus of this screening, prompting virtual screens over a large variety of head groups.

In another post, we address exactly this issue, sourcing a large number of possible head groups. We can load the primary head groups from that post.

head_pfx = "../../files/LNP_catalog/"
old_head_df = pd.read_csv(head_pfx + "frags_amines.csv.gz", compression="gzip")

mols2grid.display([Chem.MolFromSmiles(x) for x in old_head_df["SMILES"].tolist()],mol_col="mol", n_cols=7, n_rows=4)