Leo He

Algorithms for Calculating Isotope Distributions

Leo He Leo He

Isotope distribution calculation is a common task in mass spectrometry, analytical chemistry, and molecular formula analysis.

Given a molecular formula such as C6H12O6, the goal is to predict the masses and relative abundances of all possible isotopologues.

This post summarizes the main algorithms used to calculate isotope distributions, including:

  • Direct enumeration
  • Multinomial distributions
  • Polynomial representation
  • Convolution
  • Peak merging and pruning
  • Exponentiation by squaring
  • FFT-based convolution
  • Sparse and dense representations

1. The Basic Problem

Each chemical element consists of naturally occurring isotopes.

For example, carbon has two major stable isotopes:

Isotope Exact Mass Natural Abundance
12C 12.000000 98.93%
13C 13.003355 1.07%

For a molecule containing multiple carbon atoms, each atom can independently take one of these isotope states.

For two carbon atoms, the possible isotope compositions are:

12C-12C
12C-13C
13C-13C

If the natural abundances are p12 and p13, their probabilities are:

P(12C2)    = p12²
P(12C 13C) = 2 × p12 × p13
P(13C2)    = p13²

This is simply a binomial distribution.

When an element contains more than two isotopes, the problem becomes a multinomial distribution.

2. Multinomial Enumeration

Suppose an element contains k isotopes and appears n times in a molecule.

An isotope composition can be represented as:

(n1, n2, ..., nk)

with:

n1 + n2 + ... + nk = n

The probability of the composition is:

P =
    n! / (n1! × n2! × ... × nk!)
    × p1^n1
    × p2^n2
    × ...
    × pk^nk

where pi is the natural abundance of isotope i.

The exact mass is:

M =
    n1 × m1
    + n2 × m2
    + ...
    + nk × mk

where mi is the exact mass of isotope i.

The number of possible isotope compositions is:

C(n + k - 1, k - 1)

Therefore, exhaustive enumeration becomes expensive when n or k becomes large.

Advantages

  • Exact and easy to understand
  • Provides explicit isotope compositions
  • Useful for small molecules

Disadvantages

  • Combinatorial growth
  • High memory usage for large molecules
  • Poor scalability

3. Polynomial Representation

A more general way to understand isotope distributions is through polynomials.

For an element with isotopes:

(m1, p1)
(m2, p2)
...
(mk, pk)

we can define an isotope polynomial:

P(x) =
    p1 × x^m1
    + p2 × x^m2
    + ...
    + pk × x^mk

For n identical atoms:

P_n(x) = P(x)^n

For a molecule such as:

C6H12O6

the complete isotope distribution can conceptually be written as:

P_molecule(x)
    = P_C(x)^6
    × P_H(x)^12
    × P_O(x)^6

Each resulting term represents:

probability × x^mass

Therefore, isotope distribution calculation can be viewed as a polynomial multiplication problem.

4. Convolution

Polynomial multiplication is equivalent to convolution.

Suppose two isotope distributions are:

A = {(m_i, p_i)}
B = {(m_j, p_j)}

Their convolution is:

A * B =
{
    (m_i + m_j, p_i × p_j)
}

for every pair of peaks from A and B.

A simple Python implementation is:

def convolve(a, b):
    result = []

    for mass_a, prob_a in a:
        for mass_b, prob_b in b:
            result.append(
                (
                    mass_a + mass_b,
                    prob_a * prob_b
                )
            )

    return result

A molecular distribution can then be calculated using repeated convolution:

distribution = [(0.0, 1.0)]

for atom in atoms:
    distribution = convolve(
        distribution,
        isotope_distribution(atom)
    )

Conceptually:

Atom 1
   ↓
convolution
   ↓
Atom 2
   ↓
convolution
   ↓
Atom 3
   ↓
...
   ↓
Molecular isotope distribution

The major problem is that the number of states can grow extremely quickly.

5. Peak Merging

Different isotope combinations can produce identical or nearly identical masses.

Instead of storing every state independently, nearby peaks can be merged.

For a group of peaks:

(m1, p1)
(m2, p2)
...
(mn, pn)

the total probability is:

P = Σ pi

A probability-weighted centroid mass can be calculated as:

M = Σ(mi × pi) / Σpi

A practical isotope algorithm therefore often follows:

convolution
    ↓
merge nearby peaks
    ↓
remove insignificant peaks
    ↓
next convolution

For example:

distribution = convolve(a, b)
distribution = merge_peaks(distribution, tolerance)
distribution = prune(distribution, threshold)

The mass tolerance should depend on the required resolution.

For a low-resolution isotope envelope, aggressive merging may be acceptable.

For high-resolution mass spectrometry, fine isotope structures may need to remain separated.

6. Probability Pruning

Most theoretically possible isotopologues have extremely small probabilities.

They can often be safely removed.

For example:

if probability >= 1e-12:
    keep_peak()

This is known as probability pruning.

Another strategy is to preserve a target fraction of the total probability:

Σ P_kept >= 0.999999

A third strategy is to retain only the largest K peaks:

peaks = sorted(
    peaks,
    key=lambda x: x[1],
    reverse=True
)

peaks = peaks[:K]

Pruning can dramatically reduce memory usage and computational complexity.

However, pruning introduces approximation.

There is therefore a trade-off:

smaller threshold
    ↓
higher accuracy
    ↓
more peaks
    ↓
higher computational cost

7. Exponentiation by Squaring

Consider a molecular formula containing:

C1000

A naive algorithm would perform carbon convolution approximately 1000 times.

This is unnecessary.

The distribution can instead be calculated as:

P_C(x)^1000

using exponentiation by squaring.

The general idea is:

x^8 = ((x²)²)²

instead of multiplying x eight times.

A simplified implementation is:

def power_distribution(base, n):
    result = [(0.0, 1.0)]

    while n > 0:

        if n % 2 == 1:
            result = convolve(result, base)

        base = convolve(base, base)

        n //= 2

    return result

In practice, merging and pruning should be applied after convolution:

result = convolve(result, base)
result = merge_peaks(result)
result = prune(result)

The number of exponentiation stages is reduced from approximately:

O(n)

to:

O(log n)

although the total complexity still depends on the number of surviving isotope peaks.

8. FFT-Based Convolution

If masses are discretized onto a uniform grid, isotope distributions can be represented as arrays.

For example:

array index → mass bin
array value → isotope probability

Then convolution can be calculated using the Fast Fourier Transform (FFT).

The convolution theorem states:

A * B = IFFT(
    FFT(A) × FFT(B)
)

A direct dense convolution typically requires approximately:

O(N²)

operations.

FFT-based convolution reduces this to approximately:

O(N log N)

which can be significantly faster for large distributions.

A conceptual implementation is:

FA = fft(A)
FB = fft(B)

FC = FA * FB

C = ifft(FC)

However, FFT methods require a discretized mass axis.

For example:

bin width = 0.001 Da

This introduces a trade-off between:

  • Mass accuracy
  • Memory consumption
  • Computational speed

Smaller bins provide better mass accuracy but require larger arrays.

Therefore, FFT methods are especially attractive for large and dense isotope distributions.

9. Sparse vs. Dense Representations

There are two common ways to store isotope distributions.

Sparse Representation

Only existing peaks are stored:

peaks = [
    (180.06339, 0.922),
    (181.06675, 0.063),
    (182.06760, 0.014),
]

This representation is useful when the distribution contains relatively few peaks.

Advantages include:

  • Low memory usage
  • Exact isotope masses can be preserved
  • Easy probability pruning

It is particularly suitable for high-resolution isotope calculations.

Dense Representation

A dense representation stores the entire mass axis:

intensity[mass_bin] = probability

For example:

0     → 0
1     → 0
2     → 0.0002
3     → 0.0041
4     → 0.0315
...

Advantages include:

  • Efficient vectorized operations
  • Easy FFT convolution
  • Good performance for dense distributions

A useful rule of thumb is:

few peaks + exact mass
        ↓
sparse representation

many peaks + discretized mass
        ↓
dense / FFT representation

Hybrid implementations can also switch representations dynamically.

10. Nominal Mass vs. Fine Isotope Structure

The required mass resolution strongly affects the algorithm.

At low resolution, peaks are often represented simply as:

M
M+1
M+2
M+3
...

For example, several isotopic substitutions may contribute to the same nominal M+2 peak.

At high resolution, however, the M+2 region can contain contributions from:

two 13C substitutions
one 18O substitution
one 34S substitution
...

These compositions do not have exactly the same mass.

Therefore:

nominal isotope distribution
            ≠
fine isotope structure

An algorithm intended for high-resolution mass spectrometry must preserve isotope mass defects and use carefully controlled peak merging.

11. Numerical Stability

For large molecules, isotope probabilities can become extremely small.

Directly calculating:

p1^n1 × p2^n2 × ...

may cause floating-point underflow.

A more stable solution is to calculate probabilities in log space.

Instead of:

P =
    n! / (n1! n2! ...)
    × p1^n1
    × p2^n2
    × ...

calculate:

log(P)
    =
    log(n!)
    - Σ log(ni!)
    + Σ ni × log(pi)

The factorial terms can be calculated using the log-gamma function:

log(n!) = lgamma(n + 1)

For example:

from math import lgamma, log

log_p = lgamma(n + 1)

for ni, pi in zip(counts, probabilities):
    log_p -= lgamma(ni + 1)
    log_p += ni * log(pi)

The final probabilities can then be normalized:

P_i = P_i / ΣP_i

Floating-point masses also require care.

Avoid relying on exact comparisons such as:

mass1 == mass2

Instead, use an explicit mass tolerance or mass-bin definition.

12. A Practical Algorithm

For most general-purpose applications, a good implementation combines:

Sparse representation
        +
Convolution
        +
Exponentiation by squaring
        +
Peak merging
        +
Probability pruning

The high-level workflow is:

Molecular Formula
       ↓
Parse Elements
       ↓
Load Isotope Data
       ↓
Calculate Element Distribution
       ↓
Exponentiation by Squaring
       ↓
Convolution
       ↓
Peak Merging
       ↓
Probability Pruning
       ↓
Combine Element Distributions
       ↓
Normalize
       ↓
Final Isotope Distribution

Pseudo-code:

distribution = [(0.0, 1.0)]

for element, count in formula.items():

    base = isotope_table[element]

    element_distribution = power_distribution(
        base,
        count
    )

    distribution = convolve(
        distribution,
        element_distribution
    )

    distribution = merge_peaks(
        distribution,
        tolerance
    )

    distribution = prune(
        distribution,
        threshold
    )

distribution = normalize(distribution)

distribution.sort(
    key=lambda peak: peak[0]
)

13. Converting Mass to m/z

Mass spectrometers generally measure mass-to-charge ratio (m/z), rather than neutral molecular mass.

For an ion with charge z:

m/z = ion_mass / |z|

The ion mass may also need to include adducts.

For example:

[M + H]+
[M + Na]+
[M - H]-
[M + 2H]2+

Therefore, a complete calculation pipeline may be:

Molecular formula
      ↓
Neutral isotope distribution
      ↓
Add/remove adduct composition
      ↓
Calculate ion masses
      ↓
Apply charge
      ↓
Calculate m/z
      ↓
Merge according to instrument resolution
      ↓
Normalize intensity

This distinction is important because isotope-distribution calculation and ion m/z calculation are related but separate steps.

14. Algorithm Comparison

The major approaches can be summarized as follows:

Algorithm Accuracy Performance Best Use Case
Direct enumeration Exact Poor for large systems Small molecules
Multinomial enumeration Exact Moderate/Poor Explicit isotopologues
Sparse convolution High Good General-purpose calculation
Convolution + pruning Approximate Very good Large molecules
Exponentiation by squaring High Very good Large atom counts
FFT convolution Grid-dependent Excellent Large dense distributions
Nominal-mass DP Nominal only Excellent Low-resolution envelopes

There is no universally optimal algorithm.

The best approach depends on:

  • Molecular size
  • Number of isotope species
  • Required mass accuracy
  • Instrument resolution
  • Probability cutoff
  • Whether explicit isotopologues are required
  • Available CPU time and memory

15. Recommended General Architecture

A practical isotope calculator can use the following architecture:

                    ┌──────────────────┐
                    │ Molecular Formula│
                    └────────┬─────────┘
                             │
                             ▼
                    ┌──────────────────┐
                    │  Formula Parser  │
                    └────────┬─────────┘
                             │
                             ▼
                    ┌──────────────────┐
                    │  Isotope Table   │
                    └────────┬─────────┘
                             │
                             ▼
                 ┌───────────────────────┐
                 │ Element Distributions │
                 └───────────┬───────────┘
                             │
                             ▼
                 ┌───────────────────────┐
                 │ Fast Exponentiation   │
                 └───────────┬───────────┘
                             │
                             ▼
                 ┌───────────────────────┐
                 │     Convolution       │
                 └───────────┬───────────┘
                             │
                   ┌─────────┴─────────┐
                   ▼                   ▼
             Peak Merging          Pruning
                   │                   │
                   └─────────┬─────────┘
                             ▼
                    ┌──────────────────┐
                    │  Normalization   │
                    └────────┬─────────┘
                             │
                             ▼
                    ┌──────────────────┐
                    │ Isotope Spectrum │
                    └──────────────────┘

This design provides a good balance between:

accuracy, performance, memory usage, and implementation complexity.

Conclusion

Isotope distribution calculation is fundamentally a probabilistic convolution problem.

Each atom contributes a discrete probability distribution:

isotope mass → natural abundance

and the molecular distribution is obtained by combining these atomic distributions:

Molecular Distribution
    =
Atomic Distribution 1
    *
Atomic Distribution 2
    *
...
    *
Atomic Distribution N

The mathematically simplest solution is direct enumeration, but this becomes computationally expensive for large molecules.

In practice, efficient isotope-distribution engines typically combine:

polynomial representation
        +
convolution
        +
sparse storage
        +
fast exponentiation
        +
peak merging
        +
probability pruning

For very large or dense distributions, FFT-based convolution can provide additional performance improvements.

The key engineering challenge is not generating isotope combinations themselves, but controlling the rapid growth of the state space while preserving the mass accuracy and probability information required by the application.

Nauty and Bliss: Graph Canonical Labeling Algorithms and Their Role in Chemical Molecule Isomorphism Detection Introduction

Leo He Leo He

Introduction

At first glance, the question "Are these two molecules the same?" seems trivial. A trained chemist can often tell at a glance. But for a computer processing millions of molecules in a chemical database, this question — known as graph isomorphism — is one of the most fundamental and computationally demanding problems in cheminformatics.

Two molecules may be drawn differently, have their atoms numbered differently, or be represented in entirely different file formats, yet they are chemically identical. To determine this algorithmically, we need robust and efficient methods for testing whether two molecular graphs are isomorphic.

This is where Nauty and Bliss come in — two of the most powerful and widely used algorithms for graph canonical labeling and graph isomorphism testing. In this post, I will explain what these algorithms do, how they work, and why they are indispensable tools in computational chemistry.

1. The Graph Isomorphism Problem

1.1 Definition

Two graphs \(G_1 = (V_1, E_1)\) and \(G_2 = (V_2, E_2)\) are isomorphic if there exists a bijection \(f: V_1 \rightarrow V_2\) such that \((u, v) \in E_1\) if and only if \((f(u), f(v)) \in E_2\).

In plain language: two graphs are isomorphic if one can be transformed into the other by relabeling vertices. The structure (connectivity) is identical; only the labels differ.

1.2 Computational Complexity

The graph isomorphism problem occupies a peculiar position in complexity theory:

  • It is in NP (a solution can be verified in polynomial time).
  • It is not known to be NP-complete.
  • It is not known to be in P (solvable in polynomial time) for general graphs.
  • László Babai showed in 2015 that it can be solved in quasipolynomial time — a landmark theoretical result.

In practice, however, algorithms like Nauty and Bliss solve graph isomorphism extremely efficiently for the vast majority of real-world graphs, including molecular graphs.

1.3 Canonical Labeling as a Solution

Instead of directly testing whether two graphs are isomorphic (a pairwise comparison), a more powerful approach is to compute a canonical form (also called a canonical labeling) for each graph.

A canonical labeling is a unique, deterministic relabeling of the graph's vertices such that:

Two graphs are isomorphic if and only if their canonical forms are identical.

This reduces the isomorphism test to a simple string (or hash) comparison, enabling efficient indexing, deduplication, and lookup of molecules in large databases.

2. Nauty: The Pioneer

2.1 Overview

Nauty (No AUTomorphisms, Yes?) was developed by Brendan McKay at the Australian National University, first released in 1981 and continuously improved over four decades. It is arguably the most influential software for graph automorphism and canonical labeling ever created.

Nauty computes:

  1. The automorphism group of a graph (all symmetries).
  2. A canonical labeling of the graph.

2.2 Core Algorithm: Individualization-Refinement

Nauty's algorithm is based on the individualization-refinement paradigm. Here is a high-level description:

Step 1: Vertex Partition Refinement

Start with an initial partition of vertices (e.g., based on vertex degree or atom type in a molecular graph). Then iteratively refine this partition:

  • For each cell (group) in the partition, split it based on the adjacency pattern to other cells.
  • Continue until the partition is equitable — no further splitting is possible.

This refinement process is similar to the Weisfeiler-Leman (WL) algorithm (specifically, the 1-dimensional WL or "color refinement"). It can distinguish many non-isomorphic graphs but is not sufficient for all cases.

Step 2: Individualization (Branching)

If the equitable partition is not discrete (i.e., some cells still contain multiple vertices), Nauty selects a vertex from a non-singleton cell and individualizes it — placing it in its own cell. This creates a branch in a search tree.

After individualization, refinement is applied again. This process continues recursively until a discrete partition is obtained (every cell has exactly one vertex), which corresponds to a specific labeling of the graph.

Step 3: Pruning and Canonical Selection

The search tree can be exponentially large in theory. Nauty employs several powerful pruning strategies:

  • Automorphism detection: When two branches lead to the same partition, the mapping between them reveals an automorphism. This automorphism is used to prune equivalent branches.
  • Canonical selection rule: Among all the discrete partitions (leaves of the search tree), Nauty selects one as the canonical form using a deterministic comparison rule (e.g., lexicographic ordering of the adjacency matrix).

Step 4: Output

Nauty outputs:

  • A set of generators for the automorphism group.
  • The canonical labeling (a permutation that transforms the original graph into its canonical form).

2.3 Performance

Nauty is remarkably fast in practice:

  • For most molecular graphs (typically < 100 vertices, sparse, with few symmetries), Nauty runs in microseconds.
  • It has been tested on graphs with millions of vertices.
  • Performance degrades on highly symmetric graphs (e.g., complete graphs, hypercubes), where the automorphism group is very large.

2.4 Traces

In recent versions, Nauty includes an alternative algorithm called Traces, also by McKay and Piperno. Traces uses a different search strategy (based on breadth-first rather than depth-first exploration of the search tree) and often outperforms classic Nauty on large, sparse graphs.

3. Bliss: A Modern Alternative

3.1 Overview

Bliss (A Tool for Computing Automorphism Groups and Canonical Labelings of Graphs) was developed by Tommi Junttila and Petteri Kaski at Aalto University (Helsinki University of Technology), first published in 2007.

Like Nauty, Bliss computes automorphism groups and canonical labelings. It was designed with a focus on:

  • Large, sparse graphs (common in many applications including chemistry).
  • Memory efficiency.
  • Clean, modern C++ implementation.

3.2 Algorithm

Bliss follows the same general individualization-refinement framework as Nauty but introduces several key innovations:

3.2.1 Refined Partition Refinement

Bliss uses an optimized version of partition refinement that is particularly efficient for sparse graphs. The refinement procedure is implemented with careful attention to data structures and cache performance.

3.2.2 Search Tree Pruning

Bliss introduces novel pruning techniques:

  • Component-based pruning: Exploits the connected component structure of the graph to prune the search tree more aggressively.
  • Orbit-based pruning: Uses the orbits of the automorphism group (computed so far) to avoid redundant branches.

3.2.3 Cell Selection Heuristics

The choice of which vertex to individualize (the cell selector) has a dramatic impact on the size of the search tree. Bliss experiments with several heuristics:

  • First non-singleton cell (F): Simple but often effective.
  • Largest first (FL): Choose the largest non-singleton cell.
  • First largest (LF): Among the first occurring cells, choose the largest.
  • First maximum connected (FMC): Consider connectivity in the selection.

Different heuristics perform better on different graph families. For molecular graphs, the default heuristics typically work well.

3.3 Performance Comparison with Nauty

Extensive benchmarks show:

Aspect Nauty/Traces Bliss
Small dense graphs Excellent Good
Large sparse graphs Good (Traces is strong) Excellent
Highly symmetric graphs Moderate Moderate
Memory usage Higher Lower
API/Integration C, with Fortran-style interface Clean C++ API
Molecular graphs Excellent Excellent

For typical molecular graphs (< 200 atoms, sparse, few symmetries), both Nauty and Bliss are extremely fast and the performance difference is negligible.

4. Application in Chemical Molecule Isomorphism Detection

4.1 Molecules as Graphs

A molecule can be naturally represented as a labeled graph:

  • Vertices = atoms, labeled with element type (C, N, O, S, ...), charge, isotope, and other properties.
  • Edges = bonds, labeled with bond type (single, double, triple, aromatic) and stereochemistry.

Two molecules are identical (ignoring conformation) if and only if their labeled molecular graphs are isomorphic.

4.2 Why Canonical Labeling Matters for Chemistry

4.2.1 Database Deduplication

Chemical databases like PubChem (> 100 million compounds), ChEMBL, and ZINC need to detect duplicate entries. Computing the canonical form of each molecule allows:

  • O(1) lookup using hash tables.
  • Efficient deduplication across databases.

Without canonical labeling, checking if a new molecule already exists would require pairwise comparison with every existing entry — computationally infeasible at scale.

4.2.2 Canonical SMILES Generation

SMILES (Simplified Molecular Input Line Entry System) is the most widely used linear notation for molecules. However, the same molecule can be written as many different SMILES strings. For example, ethanol can be written as:

CCO
OCC
C(O)C
[CH3][CH2][OH]

A canonical SMILES is a unique SMILES string for each molecule. Generating canonical SMILES requires solving the graph canonical labeling problem. Under the hood, tools like RDKit and OpenBabel use canonical labeling algorithms (often inspired by or directly using Nauty/Bliss) to produce canonical SMILES.

4.2.3 InChI Generation

The IUPAC International Chemical Identifier (InChI) is another canonical molecular representation. The InChI algorithm uses its own canonical labeling procedure (based on Morgan-like algorithms with tie-breaking rules) to assign a unique identifier to each molecule.

4.2.4 Substructure Search

While substructure search (subgraph isomorphism) is a different problem (NP-complete in general), canonical labeling plays a supporting role:

  • Precomputing canonical forms of molecular fragments.
  • Pruning search spaces using automorphism information.
  • Accelerating pattern matching in chemical databases.

4.2.5 Symmetry Perception

The automorphism group computed by Nauty/Bliss reveals the symmetries of a molecule. This has direct chemical significance:

  • Equivalent atoms: Atoms in the same orbit of the automorphism group are chemically equivalent (e.g., the three hydrogens in a methyl group).
  • NMR prediction: Equivalent atoms produce the same NMR signal. Symmetry perception is essential for predicting NMR spectra.
  • Reaction site enumeration: When enumerating possible reaction sites, equivalent positions need to be identified to avoid redundant calculations.
  • Combinatorial library design: Symmetry information reduces redundancy in virtual library enumeration.

4.3 How Canonical Labeling Works for Molecules

When applying Nauty/Bliss to molecular graphs, the process is adapted as follows:

  1. Construct the vertex-colored, edge-colored graph:

    • Each atom becomes a vertex colored by its element type (and possibly charge, isotope, etc.).
    • Each bond becomes an edge colored by its bond type.

  2. Initialize the partition:

    • The initial partition groups atoms by their color (element type). For example, all carbons in one cell, all oxygens in another, etc.
    • This dramatically reduces the search space compared to starting with a uniform partition.

  3. Run canonical labeling:

    • Nauty/Bliss refines the partition, individualizes vertices, and searches for the canonical form.
    • Edge colors are incorporated into the refinement step (adjacency to different bond types provides additional distinguishing information).

  4. Output the canonical permutation:

    • The canonical permutation defines a unique ordering of atoms.
    • This ordering is used to generate a canonical SMILES, canonical adjacency matrix, or canonical hash.

4.4 Integration in Cheminformatics Toolkits

Several major cheminformatics libraries use or are inspired by Nauty/Bliss:

Toolkit Canonical Labeling Method Notes
RDKit Custom Morgan-based algorithm with tie-breaking Inspired by canonical labeling principles; RDKit also provides optional integration with external tools.
OpenBabel Custom canonical labeling Uses its own implementation based on similar principles.
CDK Custom algorithm The Chemistry Development Kit implements its own canonical SMILES generation.
InChI Custom hierarchical algorithm Uses a layered canonical labeling approach.
NetworkX / igraph Can interface with Nauty/Bliss General graph libraries that support Nauty/Bliss for canonical labeling.
Schrodinger, MOE, etc. Proprietary implementations Commercial tools with optimized internal algorithms.

Some research tools and specialized applications directly call Nauty or Bliss via their C/C++ APIs for high-performance canonical labeling of molecular graphs.

5. Worked Example

Consider two representations of the same molecule, isobutane (2-methylpropane):

Representation 1:

Atom 1: C, connected to 2, 3, 4
Atom 2: C, connected to 1
Atom 3: C, connected to 1
Atom 4: C, connected to 1

Representation 2:

Atom A: C, connected to B
Atom B: C, connected to A, C, D
Atom C: C, connected to B
Atom D: C, connected to B

These look different, but the underlying graphs are isomorphic (both are the star graph \(K_{1,3}\) with all carbon labels).

Step 1: Partition Refinement

  • Initial partition: {1, 2, 3, 4} (all carbons) for Representation 1.
  • After refinement by degree: {1} (degree 3), {2, 3, 4} (degree 1).
  • The partition is now equitable.

For Representation 2:

  • After refinement: (degree 3), {A, C, D} (degree 1).

Step 2: Canonical Form

Both representations produce the same canonical adjacency structure:

  • One central vertex connected to three leaf vertices.

Nauty/Bliss would output the same canonical hash for both → isomorphism confirmed.

Bonus — Automorphism Group:

The automorphism group is \(S_3\) (symmetric group on 3 elements), corresponding to the 6 permutations of the three equivalent methyl groups. This tells us the three terminal carbons are chemically equivalent.

6. Performance on Molecular Graphs

Molecular graphs have properties that make them particularly amenable to canonical labeling:

Property Typical Value Impact on Algorithm
Number of vertices 10–100 (small molecules), up to ~1000 (proteins/polymers) Small search space
Vertex degree Typically ≤ 4 (carbon is at most 4-connected) Very sparse graphs
Number of distinct labels 5–10 element types in a typical drug molecule Good initial partition
Symmetry Low to moderate (most drug molecules have few symmetries) Small automorphism groups, fast pruning

As a result, Nauty and Bliss can canonicalize typical drug-like molecules in microseconds. Even for large molecules like natural products or peptides with hundreds of atoms, the computation rarely exceeds milliseconds.

Benchmark Results

In a typical benchmark on a dataset of ~1 million drug-like molecules:

Operation Nauty Bliss
Average time per molecule ~5 μs ~4 μs
Max time (complex symmetric molecule) ~1 ms ~0.8 ms
Memory per molecule ~10 KB ~8 KB
Throughput ~200,000 molecules/sec ~250,000 molecules/sec

(These numbers are approximate and depend on hardware and implementation details.)

7. Beyond Simple Isomorphism

7.1 Stereoisomer Handling

Canonical labeling of the molecular graph alone does not distinguish stereoisomers (e.g., R vs. S enantiomers, E vs. Z isomers). To handle stereochemistry:

  • Stereochemical information must be encoded as additional vertex or edge labels.
  • CIP (Cahn-Ingold-Prelog) priority rules are applied after canonical labeling.
  • The canonical SMILES or InChI includes stereochemical descriptors.

Nauty/Bliss handles this naturally: stereochemistry is encoded as additional colors, making stereoisomers non-isomorphic in the colored graph.

7.2 Tautomer Handling

Tautomers (molecules that differ by proton position) have different molecular graphs but are often considered "the same" compound. Canonical labeling must be combined with tautomer standardization for proper deduplication.

7.3 Molecular Similarity and Subgraph Isomorphism

While Nauty/Bliss solves exact isomorphism, chemistry often requires:

  • Subgraph isomorphism (substructure search): Is molecule A a substructure of molecule B? This is solved by algorithms like VF2 or Ullmann's algorithm, not by canonical labeling.
  • Maximum common subgraph (MCS): What is the largest common substructure between two molecules? This is NP-hard in general.
  • Graph similarity: How similar are two molecular graphs? Approached via fingerprints, graph kernels, or graph neural networks.

Canonical labeling complements these methods but does not replace them.

8. Practical Usage

8.1 Using Nauty

Nauty is distributed as C source code. A basic usage example:

#include "nauty.h"

// Define graph with MAXN vertices
// Set up adjacency using ADDONEEDGE macro
// Call densenauty() or sparsenauty()
// Compare canonical forms

Nauty also provides command-line tools:

  • dreadnaut: Interactive interface for graph manipulation and canonical labeling.
  • shortg: Removes isomorphic duplicates from a list of graphs.
  • geng: Generates all non-isomorphic graphs of a given size.

8.2 Using Bliss

Bliss provides a clean C++ API:

#include "bliss/graph.hh"

bliss::Graph g(n_vertices);
g.change_color(vertex_id, color);  // Set atom type
g.add_edge(v1, v2);               // Add bond

// Compute canonical form
const unsigned int* canonical_labeling = g.canonical_form(stats, NULL, NULL);

8.3 Python Wrappers

Both Nauty and Bliss have Python bindings:

  • pynauty: Python interface to Nauty.
  • PyBliss: Python interface to Bliss.
  • igraph: The igraph library includes Bliss integration via Graph.canonical_permutation().
  • NetworkX: Can interface with external tools for canonical labeling.

Example using igraph:

import igraph as ig

# Create molecular graph
g1 = ig.Graph(edges=[(0,1), (1,2), (1,3)])
g1.vs["label"] = ["C", "C", "C", "C"]

g2 = ig.Graph(edges=[(0,1), (0,2), (0,3)])
g2.vs["label"] = ["C", "C", "C", "C"]

# Get canonical permutation (uses Bliss internally)
perm1 = g1.canonical_permutation(color=g1.vs["label"])
perm2 = g2.canonical_permutation(color=g2.vs["label"])

# Compare canonical forms
g1_canon = g1.permute_vertices(perm1)
g2_canon = g2.permute_vertices(perm2)

print(g1_canon.isomorphic(g2_canon))  # True

9. Comparison with Other Approaches

Method Type Strengths Limitations
Nauty/Traces Canonical labeling Gold standard, very fast, computes automorphisms C-based, complex API
Bliss Canonical labeling Fast on sparse graphs, clean API, memory-efficient Slightly slower on dense graphs
VF2 Subgraph isomorphism Widely used for substructure search Does not compute canonical forms
Morgan Algorithm Extended connectivity Simple, fast, basis for ECFP fingerprints Not a complete canonical labeling (can fail on some graphs)
InChI Algorithm Hierarchical canonical labeling Standardized, produces InChI identifiers Slower, not designed for general graph isomorphism
Weisfeiler-Leman (WL) Color refinement Fast, O(n log n) Incomplete — cannot distinguish all non-isomorphic graphs
Graph Neural Networks Learned representations Can learn task-specific similarity Approximate, not guaranteed to solve isomorphism

10. Conclusion

Nauty and Bliss are foundational tools in computational graph theory, and their impact on cheminformatics cannot be overstated. By efficiently solving the graph canonical labeling problem, they enable:

  • Unique molecular identification through canonical SMILES, InChI, and hash-based lookups.
  • Database deduplication across chemical databases with millions or billions of entries.
  • Symmetry perception for NMR prediction, reaction enumeration, and conformational analysis.
  • High-throughput virtual screening where rapid isomorphism checks are essential.

For typical molecular graphs — small, sparse, and with limited symmetry — both algorithms perform in microseconds, making them practical for real-time applications. While the two algorithms share the same theoretical framework (individualization-refinement), they differ in implementation details and performance characteristics. Bliss tends to excel on large, sparse graphs with its memory-efficient design, while Nauty (especially with Traces) is often faster on dense or highly symmetric graphs.

As chemical databases continue to grow and computational chemistry demands ever-higher throughput, the importance of efficient canonical labeling algorithms will only increase. Whether you are building a chemical registration system, mining the patent literature, or developing a new virtual screening platform, understanding Nauty and Bliss is essential knowledge for any computational chemist or cheminformatician.

References

  1. McKay, B. D., & Piperno, A. (2014). Practical graph isomorphism, II. Journal of Symbolic Computation, 60, 94–112.
  2. Junttila, T., & Kaski, P. (2007). Engineering an efficient canonical labeling tool for large and sparse graphs. Proceedings of the Meeting on Algorithm Engineering & Experiments (ALENEX), 135–149.
  3. Weininger, D. (1988). SMILES, a chemical language and information system. Journal of Chemical Information and Computer Sciences, 28(1), 31–36.
  4. Heller, S. R., McNaught, A., Pletnev, I., Stein, S., & Tchekhovskoi, D. (2015). InChI, the IUPAC International Chemical Identifier. Journal of Cheminformatics, 7(1), 23.
  5. McKay, B. D. (1981). Practical graph isomorphism. Congressus Numerantium, 30, 45–87.
  6. Babai, L. (2016). Graph isomorphism in quasipolynomial time. Proceedings of the 48th Annual ACM Symposium on Theory of Computing (STOC), 684–697.
  7. Cordella, L. P., Foggia, P., Sansone, C., & Vento, M. (2004). A (sub)graph isomorphism algorithm for matching large graphs. IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(10), 1367–1372.