Merck KGaA: Quantum-Classical Hybrid Workflows for Fragment-Based Drug Design

The Problem

Fragment-based drug design (FBDD) is a widely used approach in pharmaceutical research. Instead of searching a library of complete drug-sized molecules, medicinal chemists screen thousands of small molecular fragments that each make simple, weak contact with a target protein. Promising fragments are then combined and grown into larger molecules with potent binding.

The combinatorial challenge: given a set of 500 to 2000 fragment candidates, which combinations produce the best drug lead? Two sub-problems arise. First, each fragment’s binding affinity requires an energy calculation. Classical methods (molecular mechanics, MM-GBSA) are fast but approximate; quantum chemistry methods (DFT, CCSD) are accurate but expensive. Second, assembling fragments requires solving a combinatorial selection problem subject to chemical compatibility constraints.

Merck KGaA’s quantum computing lab, based in Darmstadt, addressed both sub-problems with a hybrid workflow: VQE for binding energies and QAOA for the combinatorial assembly step.

QAOA for Fragment Selection

The fragment selection problem can be framed as a graph problem. Fragments are nodes. Edges carry two types of weights: positive synergy (fragments that cooperate in binding) and negative penalty (steric clashes or chemical incompatibility). The objective is to select a subset of fragments that maximizes total synergy while respecting compatibility constraints.

from qiskit_optimization import QuadraticProgram
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler
import numpy as np

# Fragment selection as a quadratic program
# Variables: x_i = 1 if fragment i is selected
n_fragments = 12

# Simulated synergy scores between fragment pairs (upper triangle)
np.random.seed(42)
synergy = np.random.uniform(-2, 5, (n_fragments, n_fragments))
synergy = np.triu(synergy, k=1)

# Build quadratic program
qp = QuadraticProgram(name="fragment_selection")
for i in range(n_fragments):
    qp.binary_var(name=f"x{i}")

# Objective: maximize synergy (minimize negative synergy)
linear_terms = {}
quadratic_terms = {}

for i in range(n_fragments):
    for j in range(i + 1, n_fragments):
        if abs(synergy[i, j]) > 0.1:
            quadratic_terms[(f"x{i}", f"x{j}")] = -synergy[i, j]

# Constraint: select at most 4 fragments (budget constraint)
qp.minimize(linear=linear_terms, quadratic=quadratic_terms)
qp.linear_constraint(
    linear={f"x{i}": 1 for i in range(n_fragments)},
    sense="<=",
    rhs=4,
    name="fragment_budget"
)

print(qp.export_as_lp_string())

# Solve with QAOA
sampler = Sampler()
qaoa = QAOA(sampler=sampler,
            optimizer=COBYLA(maxiter=300),
            reps=2)

optimizer = MinimumEigenOptimizer(qaoa)
result = optimizer.solve(qp)

print(f"Selected fragments: {[i for i in range(n_fragments) if result.x[i] > 0.5]}")
print(f"Objective value: {result.fval:.3f}")

VQE for Fragment Binding Energy

The binding energy component uses Qiskit Nature to set up electronic structure Hamiltonians for small fragment-protein subsystems. For fragments with 3-5 heavy atoms, the active space can be reduced to 4-8 qubits using a minimal basis set and active space selection.

The workflow: (1) generate a classical geometry for the fragment in the protein binding site using molecular dynamics, (2) extract an active space Hamiltonian using PySCF via Qiskit Nature, (3) run VQE with a UCCSD-inspired ansatz, (4) return the energy estimate to the fragment selection scoring function.

For the fragments tested (acetamide, imidazole, and pyrazole derivatives binding a serine protease), VQE energies matched MM-GBSA estimates within 0.5 kcal/mol for relative binding rankings, which is the operationally relevant metric for fragment prioritization.

Results

Merck benchmarked the hybrid workflow against fully classical pipelines on a curated set of 20 fragments targeting a kinase active site:

• QAOA fragment selection (12 fragments, p=2): Matched classical genetic algorithm ranking on 8 of 10 test cases; 2 cases had equivalent fragment combinations ranked differently with similar total synergy scores
• VQE binding energies (4-8 qubit active spaces): Correlation coefficient r=0.89 against classical MM-GBSA on relative binding free energies across tested fragments
• Runtime: VQE on simulator was 10-50x slower than MM-GBSA; QAOA classical overhead (kernel matrix, optimization) dominated over quantum circuit execution time at this scale
• Scalability gap: Drug-relevant fragments (10-20 heavy atoms) require 20-40 active space orbitals, mapping to 40-80 qubits. Fault-tolerant hardware with 100+ logical qubits is the estimated threshold for genuine advantage

What’s Next

Merck KGaA published the methodology in a white paper and submitted experimental results to a peer-reviewed journal. The roadmap has three phases:

1. Near term: Expand VQE experiments to fragments with up to 8 heavy atoms using IBM Quantum’s 127-qubit Eagle processors with error mitigation
2. Mid term: Test full fragment-to-lead workflows on early fault-tolerant devices as they become available through IBM’s roadmap
3. Long term: Target protein-ligand binding calculations at drug-relevant scales (100+ logical qubits), aiming to replace or validate classical scoring functions in prospective drug discovery campaigns

Learn more: Qiskit Reference