Drug Discovery with Machine Learning

Suppose a pharmaceutical research group wants to find a molecule that binds tightly to a specific viral protein — say, a protease that a virus needs to cut its own proteins into working pieces. The chemists have access to a well-stocked compound library and a robotic high-throughput screening (HTS) rig that can run roughly 100,000 binding assays a day, which is a genuinely fast, realistic industrial rate. The catch is the size of the space they are searching. Chemists estimate that the universe of "drug-like" small molecules — organic compounds under about 500 daltons that could plausibly be turned into a pill — numbers somewhere around 10⁶⁰, an estimate first popularised by Bohacek, McMartin and Guida in 1996 and still the order-of-magnitude figure the field cites. At 100,000 assays a day, clearing that library by brute force would take roughly 10⁵⁵ days, or about 2.7 × 10⁵² years. The universe is about 1.38 × 10¹⁰ years old. Exhaustive wet-lab screening of chemical space is not merely expensive — it is off by roughly 42 orders of magnitude from feasible, even before you account for the fact that a compound also has to be synthesisable, non-toxic, and orally absorbable. This is the actual problem machine learning is hired to solve in drug discovery: not "invent chemistry," but replace blind enumeration with a learned function that ranks candidates so the lab only ever tests the few thousand most promising ones.

That reframing is what this chapter teaches. You already know how to represent a network as a graph and traverse it; you already know regression, feature vectors, and how a neural network layer transforms one representation into another. Drug discovery with ML is, mechanically, an application of exactly those tools to a domain where the "graph" is a molecule and the "label" is a measured biological effect. The rest of this chapter builds that mapping precisely, works a full numeric example by hand, and is honest about where the technique stops and the biology begins.

A molecule is a graph you already know how to process

Every organic molecule is, structurally, exactly the kind of object you studied in your data structures course: a set of vertices (atoms) connected by a set of edges (bonds), each edge carrying a weight-like attribute (bond order — single, double, aromatic). Chemists have a compact linear encoding for this graph called SMILES (Simplified Molecular Input Line Entry System). Aspirin, for instance, is written as CC(=O)Oc1ccccc1C(=O)O — the letters are atoms, parentheses mark branches, digits mark ring closures, and lowercase letters mark aromatic atoms. A SMILES parser turns that string into an adjacency structure in the same way you would parse a bracket expression into a tree: a stack tracks open branches and ring-closure digits are matched like parentheses to add the missing edge. Once parsed, the molecule is a small, sparse, undirected, labelled graph — typically 15 to 60 atoms for a drug-like compound, with each atom carrying a feature vector (element, formal charge, hybridisation, aromaticity, number of attached hydrogens) and each bond carrying its order.

This matters pedagogically because it means every graph algorithm you already know — breadth-first search, connected components, shortest path — applies unmodified to a molecule. A ring-detection routine is just cycle detection. Finding the longest conjugated chain is a longest-path computation on a subgraph. And, as you'll see in the next two sections, the two dominant ML approaches to molecular property prediction are best understood as, respectively, "flatten the graph into a fixed-length vector and run classical regression on it" (fingerprints and QSAR) and "run a BFS-like traversal where each node accumulates a real-valued vector instead of a visited flag" (graph neural networks).

Flattening the graph: fingerprints and QSAR

The oldest ML approach to drug discovery, still heavily used, sidesteps the graph structure by hashing it into a fixed-length bit vector called a molecular fingerprint. The most common family, circular fingerprints (Extended-Connectivity Fingerprints, or ECFP, introduced by Rogers and Hahn in 2010), work by a procedure that is, again, a graph algorithm you know: for each atom, hash its own identity into an integer; then, for a chosen radius r, repeatedly hash each atom's current code together with the sorted codes of its neighbours to get a new code — exactly one round of a message-passing update, except the "aggregation function" is a hash instead of a sum, and there is no learning involved. Every hashed code that appears, across every atom and every radius up to r, is folded (via modulo) into a bit position in a fixed-length vector, typically 1024 or 2048 bits long. Two molecules that share substructures light up overlapping bits; the Tanimoto (Jaccard) similarity between their bit vectors becomes a cheap structural-similarity score.

Quantitative Structure-Activity Relationship (QSAR) modelling then treats that fingerprint — or a smaller, hand-picked set of physicochemical descriptors such as molecular weight, calculated octanol-water partition coefficient (logP), hydrogen-bond donor and acceptor counts — as an ordinary feature vector x, and regresses it against a measured biological activity y (commonly pIC50, the negative log of the concentration needed to inhibit a target by 50%) using linear regression, random forests, or gradient boosting. This is precisely the regression pipeline from your ML foundations, with molecules as rows. One famous descriptor-based rule that predates ML but shows the same logic is Lipinski's Rule of Five (1997): oral drugs tend to satisfy molecular weight ≤ 500, logP ≤ 5, hydrogen-bond donors ≤ 5, and acceptors ≤ 10 — a crude but still-used linear filter on exactly this kind of feature vector.

QSAR's limitation is that a fixed hand-designed fingerprint throws away information the model might have used — two structurally different scaffolds that happen to hash into similar bit patterns are treated as similar, and vice versa. Graph neural networks exist to fix this by learning the feature-extraction step end to end instead of relying on a fixed hash.

Graph neural networks: message passing is BFS with vectors instead of flags

Recall breadth-first search: you initialise every node's state to "unvisited," then repeatedly look at each node's neighbours and update state based on what has reached it so far. A message-passing graph neural network (the architecture underlying nearly every modern molecular-property model, from Google's message-passing neural networks for chemistry to the graph models inside Insilico Medicine's and Iktos's drug-design pipelines) runs the identical control structure, but instead of a boolean flag, every node carries a real-valued embedding vector, and instead of "mark visited," each round applies a differentiable update.

Formally, at layer k, every atom v holds an embedding h_v^(k). One message-passing layer computes:

h_v^(k+1) = ReLU( W^(k) · aggregate( h_v^(k), {h_u^(k) : u ∈ N(v)} ) )

where N(v) is the set of atoms bonded to v — exactly the adjacency-list neighbours a BFS would enqueue. A common, simple choice of aggregate (used in the Graph Convolutional Network of Kipf and Welling, 2017) is the mean of the node's own embedding and its neighbours' embeddings, normalised by degree. W^(k) is a learned weight matrix, shared across every atom in the molecule and every molecule in the dataset — the same parameters do the work everywhere, just as a single recursive function definition works at every node of a tree. Stack k such layers and, exactly as in BFS, atom v's final embedding depends on everything within k hops of it — its k-hop receptive field. After 3–4 layers (the typical depth used in production molecular GNNs), every atom's vector encodes information about its local chemical environment: what it's bonded to, what that's bonded to, and so on. A final "readout" step pools all the atom embeddings — commonly by mean or sum — into one fixed-length vector for the whole molecule, which a small feed-forward head then maps to a predicted property: binding affinity, solubility, toxicity, whatever the model was trained on.

Worked example: one message-passing layer, by hand

Take a deliberately small toy molecule — the three-carbon skeleton of propene, CH₂=CH–CH₃ — reduced to a path graph on three atoms C1–C2–C3, with a double bond C1=C2 and a single bond C2–C3. Give each atom a 2-dimensional initial feature vector encoding (is-sp², scaled hydrogen count):

h1(0) = [1, 2]   # =CH2, sp2, 2 attached H
h2(0) = [1, 1]   # =CH-, sp2, 1 attached H
h3(0) = [0, 3]   # -CH3, sp3, 3 attached H

Degrees in this path graph: deg(C1) = 1, deg(C2) = 2, deg(C3) = 1. Using the mean-aggregation rule — self plus neighbours, divided by (1 + degree) — the pre-activation aggregate at C2 is:

agg(C2) = (h2(0) + h1(0) + h3(0)) / (1 + deg(C2))
        = ([1,1] + [1,2] + [0,3]) / 3
        = [2, 6] / 3
        = [0.667, 2.000]

Now apply a learned linear transform. Take a small illustrative weight matrix W = [[1,-1],[0,1]] — chosen only so the arithmetic is easy to check by hand, not a trained value. Applying z = W·h to a 2-vector [a,b] gives [a−b, b], so:

z(C2) = [0.667 − 2.000, 2.000] = [-1.333, 2.000]
h2(1) = ReLU(z(C2)) = [0, 2.000]

The identical rule, applied at C1 (neighbours = {C2} only, degree 1) and C3 (neighbours = {C2} only, degree 1), gives:

agg(C1) = ([1,2]+[1,1]) / 2 = [1.000, 1.500] → z = [-0.500, 1.500] → h1(1) = [0, 1.500]
agg(C3) = ([0,3]+[1,1]) / 2 = [0.500, 2.000] → z = [-1.500, 2.000] → h3(1) = [0, 2.000]

The full computation, verified in code rather than only by hand, is:

import numpy as np

H0 = np.array([[1, 2], [1, 1], [0, 3]], dtype=float)   # h1,h2,h3 initial features
A  = np.array([[0, 1, 0],
               [1, 0, 1],
               [0, 1, 0]], dtype=float)                 # adjacency, path graph

degree = A.sum(axis=1)                                  # [1, 2, 1]
denom  = (1 + degree).reshape(-1, 1)                     # [[2],[3],[2]]
H_pre  = (H0 + A @ H0) / denom                           # aggregate self + neighbours

W = np.array([[1, -1], [0, 1]], dtype=float)
H1 = np.maximum(H_pre @ W.T, 0)                          # linear transform + ReLU

readout = H1.mean(axis=0)                                # graph-level embedding

w_out, b_out = np.array([0.5, 1.2]), -1.0                # illustrative output head
logit = readout @ w_out + b_out
score = 1 / (1 + np.exp(-logit))

print(H_pre)      # [[1.    1.5 ] [0.667 2.   ] [0.5   2.  ]]
print(H1)         # [[0.  1.5] [0. 2. ] [0.  2. ]]
print(readout)    # [0.    1.833]
print(score)      # 0.7685...

Mean-pooling the three layer-1 embeddings gives the graph-level vector [0, 1.833] — check: (1.5 + 2.0 + 2.0)/3 = 5.5/3 = 1.8333, and the first coordinate is (0+0+0)/3 = 0. Feeding that through an illustrative output head w_out = [0.5, 1.2], b = −1.0 gives a logit of 0.5·0 + 1.2·1.8333 − 1.0 = 2.2 − 1.0 = 1.2 exactly (since 1.8333… × 1.2 = 2.2 precisely, because 5.5/3 × 1.2 = 6.6/3 = 2.2). The sigmoid of 1.2 is 1/(1+e^−1.2) ≈ 1/1.3012 ≈ 0.769 — an illustrative "likely binder" score. This end-to-end chain — graph, aggregation, linear transform, ReLU, readout, output head — is the entire mechanism of a molecular GNN; a production model just uses 30–300 dimensional embeddings, 3–6 layers, and weights learned by gradient descent on millions of labelled molecules from databases like ChEMBL (a public bioactivity database) rather than the two illustrative matrices used here.

The diagram: message passing traced on this exact molecule

Where this fits in a real discovery pipeline

A trained GNN or QSAR model is a filter, not a finish line. The real pipeline runs: target identification (find the protein whose function, if blocked or activated, treats the disease) → virtual screening or generative design (rank or generate candidate molecules against that target using models like the one above) → hit-to-lead optimisation (chemists modify the top-ranked scaffolds to improve potency) → ADMET prediction (separate ML models estimate Absorption, Distribution, Metabolism, Excretion and Toxicity, since a molecule that binds perfectly but is destroyed by the liver in ten seconds is useless) → wet-lab synthesis and assay of the survivors → animal studies → human Phase 1/2/3 clinical trials. ML today accelerates the first three stages substantially and has a growing role in ADMET prediction; it has essentially no ability to shorten the trial stages, because those are rate-limited by human biology and regulatory safety requirements, not by search.

Three real cases show both the power and the limits. DeepMind's AlphaFold2, which in 2020 solved protein structure prediction to near-experimental accuracy at the CASP14 competition, now supplies predicted 3D structures for hundreds of millions of proteins through the AlphaFold Database — feeding directly into structure-based virtual screening, since docking a candidate molecule against a target requires knowing the target's 3D shape. BenevolentAI used a knowledge-graph reasoning system (a different ML technique from the GNN above — it mines relationships across biomedical literature rather than molecular graphs) to identify baricitinib, an existing rheumatoid-arthritis drug, as a plausible COVID-19 treatment by connecting its known mechanism to viral cell entry and cytokine-storm pathways; it received FDA emergency authorisation in November 2020 — a repurposing win, not a from-scratch discovery. Insilico Medicine's PandaOmics and Chemistry42 platforms generated INS018_055, a novel small molecule for idiopathic pulmonary fibrosis, and are commonly credited with taking roughly 18 months from target identification to a nominated preclinical candidate — against a historical average closer to four or five years for that stage — making it one of the first largely AI-generated molecules to reach human clinical trials. Even in that case, the compound still had to go through the full multi-year sequence of Phase 1 safety and Phase 2/3 efficacy trials in humans; the acceleration was confined to the computational discovery stage.

The misconception worth correcting

The natural misreading of everything above is: "if a model can generate and rank candidate molecules in silico, the years-long drug pipeline effectively collapses." It doesn't, for two separable reasons — one about biology, one about the model itself. The biological reason is the one already stated: clinical trials measure how a molecule behaves inside a living human immune and metabolic system over months to years, and no computation shortens that measurement, only the search that produces the candidate worth measuring. The modelling reason is a real weakness of the fingerprint- and graph-based similarity that underlies these systems: the activity cliff. Two molecules that differ by a single atom or a single stereo-centre can have wildly different biological activity — one a potent inhibitor, the near-identical twin inert or even toxic — while their fingerprints or GNN embeddings, which are built from local structural similarity, place them almost on top of each other. A model trained mostly on one chemical series will also generalise poorly to a genuinely novel scaffold outside its training distribution, precisely the "out-of-distribution" failure mode from your general ML foundations, here with the added twist that the training data (measured bioactivities in ChEMBL and similar databases) is itself skewed toward chemical families that medicinal chemists have historically explored, not toward the full 10⁶⁰-molecule space. The honest summary: ML compresses the hit-finding stage from years to months and helps triage which candidates deserve a chemist's and a clinician's time; it does not remove the need for wet-lab validation, and most computationally flagged candidates still fail downstream, exactly as most candidates always have.

Active recall

Attempt each question before reading its answer.

1. A biotech's HTS robot can test 200,000 compounds a day. Using the ~10⁶⁰ drug-like chemical-space estimate, roughly how many years would exhaustive screening take, and how does that compare to the age of the universe (≈1.38×10¹⁰ years)?
2. In the SMILES string CC(=O)Oc1ccccc1C(=O)O (aspirin), what do the digit "1"s and the parentheses each represent structurally?
3. A water-like star graph has a central atom O bonded to two atoms H_a and H_b (O has degree 2; each H has degree 1). Given h_O⁽⁰⁾=[3.5, 2], h_Ha⁽⁰⁾=[2.1, 0], h_Hb⁽⁰⁾=[2.1, 0], and W=[[2,0],[1,1]], compute h_O⁽¹⁾ using the same mean-aggregation-then-ReLU rule used for C2 in the worked example.
4. Why does a high Tanimoto fingerprint similarity between two candidate molecules NOT guarantee similar measured biological activity? Name the phenomenon.
5. Using the readout vector [0, 1.833] from the worked example, w_out=[0.5, 1.2], and b=−1.0, verify the final predicted score is ≈0.77, then explain in one sentence why a high score here does not mean the molecule is ready to prescribe.
6. QSAR regression and a graph neural network both eventually produce a fixed-length vector representation of a molecule before predicting activity. What is the essential difference in how that vector is produced?

Answers.

1. 10⁶⁰ ÷ 2×10⁵ = 5×10⁵⁴ days ÷ 365.25 ≈ 1.37×10⁵² years — about 10⁴² times the age of the universe. Brute-force screening is not merely slow, it is physically impossible within any realistic timeframe, which is the actual justification for using a learned ranking function instead of exhaustive testing.

2. The "1"s are ring-closure labels: the first atom marked "1" and the later atom marked "1" are bonded to each other, closing the benzene ring, in exactly the way a parser matches a pair of bracket tokens. The parentheses mark a branch off the main chain — here, (=O) attaches a double-bonded oxygen to the preceding carbon, then parsing resumes from that same carbon once the branch closes.

3. Aggregate: (h_O+h_Ha+h_Hb)/(1+deg(O)) = ([3.5+2.1+2.1, 2+0+0])/3 = [7.7, 2.0]/3 = [2.567, 0.667]. Transform with W=[[2,0],[1,1]]: z = [2×2.567 + 0×0.667, 1×2.567 + 1×0.667] = [5.133, 3.233]. Both entries are already positive, so ReLU leaves them unchanged: h_O⁽¹⁾ = [5.133, 3.233].

4. Activity cliffs. Fingerprint similarity is a measure of local structural overlap, but biological activity depends on precise 3D shape and electronic complementarity with a binding pocket; a single substituent swap can flip a key steric or electronic interaction while barely moving the fingerprint, so structural closeness is not a reliable proxy for activity closeness.

5. logit = 0.5×0 + 1.2×1.833 − 1.0 = 2.2 − 1.0 = 1.2; σ(1.2) = 1/(1+e^−1.2) ≈ 1/1.301 ≈ 0.769, matching the diagram. A high predicted-binding score only says the model expects the molecule to interact with the target; it says nothing about toxicity, metabolic stability, or how the molecule behaves in a human body, all of which are only established through the ADMET, animal, and clinical-trial stages that follow.

6. QSAR uses a fixed, hand-designed procedure (a hashing scheme for fingerprints, or manually chosen descriptors like molecular weight and logP) to produce the vector before any learning happens — only the final regression step is trained. A GNN learns the vector itself: the weight matrices used in every message-passing round are trained end to end from labelled data, so the network can discover which structural patterns matter for the specific property being predicted, rather than relying on a hash function designed without that property in mind.

Think About It

Think about this: How would you explain drug discovery with machine learning to a friend who has never seen a computer? What real-world analogy would you use? Imagine you had to build a system using these concepts — what would be your first step? Try this: before moving on, write down three things you learned and one question you still have.

Practice Exercises

Now it is time to practice! Complete these challenges to solidify your understanding:

• Exercise 1: Write a short program that demonstrates the core concept from this chapter. Test it with at least 3 different inputs.
• Exercise 2: Find a real-world example where drug discovery with machine learning is used in an Indian company (like TCS, Infosys, Flipkart, or ISRO). Write a paragraph explaining the connection.
• Exercise 3: Create a mind-map connecting drug discovery with machine learning to at least 3 other topics you have studied.