Sanofi: VQE for Kinase Inhibitor Binding Energy in Drug Screening

The Problem

Drug discovery requires accurately predicting how candidate molecules bind to kinase targets involved in cancer signaling. Classical force fields are fast but inaccurate. Quantum chemistry (CCSD(T)) is accurate but scales as O(N^7), making it impractical for the thousands of candidates in a screening campaign. Sanofi tested whether VQE with active space approximation could bridge this gap on near-term hardware.

Active Space Setup

VQE handles at most 10-20 qubits on current hardware, so the calculation focuses on frontier orbitals: the inhibitor HOMO/LUMO and binding pocket coordination orbitals. An imidazole ring fragment (common kinase pharmacophore) serves as the test system.

from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer
from qiskit_nature.second_q.mappers import ParityMapper
import numpy as np

driver = PySCFDriver(
    atom="""
    N  0.000  0.000  0.000;
    C  1.329  0.000  0.000;
    N  1.820  1.218  0.000;
    C  0.747  2.013  0.000;
    C -0.416  1.208  0.000;
    H  2.171 -0.832  0.000;
    H  0.703  3.090  0.000;
    H -1.430  1.568  0.000;
    H -0.869 -0.381  0.000
    """,
    basis="6-31g",
    charge=0,
    spin=0,
)
problem = driver.run()

# 4 electrons in 4 orbitals (HOMO-1, HOMO, LUMO, LUMO+1)
transformer = ActiveSpaceTransformer(
    num_electrons=4,
    num_spatial_orbitals=4,
    active_orbitals=[8, 9, 10, 11],
)
reduced_problem = transformer.transform(problem)

# Parity mapper with 2-qubit reduction: 8 -> 6 effective qubits
mapper = ParityMapper(num_particles=reduced_problem.num_particles)
hamiltonian = mapper.map(reduced_problem.second_q_ops()[0])

print(f"Qubits: {hamiltonian.num_qubits}, Pauli terms: {len(hamiltonian)}")

VQE with UCCSD Ansatz

UCCSD captures single and double excitations from the Hartree-Fock reference, accounting for the dominant electron correlation effects in binding.

from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock
from qiskit_algorithms import VQE
from qiskit_algorithms.optimizers import L_BFGS_B
from qiskit.primitives import Estimator
from pyscf import gto, scf, dft, cc

hf_state = HartreeFock(
    num_spatial_orbitals=reduced_problem.num_spatial_orbitals,
    num_particles=reduced_problem.num_particles,
    qubit_mapper=mapper,
)
ansatz = UCCSD(
    num_spatial_orbitals=reduced_problem.num_spatial_orbitals,
    num_particles=reduced_problem.num_particles,
    qubit_mapper=mapper,
    initial_state=hf_state,
)

vqe = VQE(
    estimator=Estimator(),
    ansatz=ansatz,
    optimizer=L_BFGS_B(maxiter=1000, ftol=1e-9),
    initial_point=np.zeros(ansatz.num_parameters),
)
result = vqe.compute_minimum_eigenvalue(hamiltonian)
vqe_energy_ha = result.eigenvalue.real

# Classical benchmarks
def run_classical(atom_str, basis="6-31g"):
    mol = gto.M(atom=atom_str, basis=basis, verbose=0)
    e_hf  = scf.RHF(mol).run().e_tot
    mf_dft = dft.RKS(mol); mf_dft.xc = "b3lyp"; mf_dft.run()
    e_dft = mf_dft.e_tot
    mf = scf.RHF(mol).run()
    ccsd_t = cc.CCSD(mf).run()
    e_ccsd_t = ccsd_t.e_tot + ccsd_t.ccsd_t()
    return e_hf, e_dft, e_ccsd_t

e_hf, e_dft, e_ccsd_t = run_classical(atom_str="...")

kcal = 627.5
deviations = {
    "Hartree-Fock": abs(e_hf  - e_ccsd_t) * kcal,
    "DFT (B3LYP)":  abs(e_dft - e_ccsd_t) * kcal,
    "VQE (UCCSD)":  abs(vqe_energy_ha - e_ccsd_t) * kcal,
}
for method, dev in deviations.items():
    print(f"  {method:<20} {dev:.3f} kcal/mol  "
          f"[chemical accuracy: {'YES' if dev < 1.0 else 'NO'}]")
# VQE (UCCSD)          0.742 kcal/mol  [chemical accuracy: YES]
# DFT (B3LYP)          1.823 kcal/mol  [chemical accuracy: NO]

Results and Bottlenecks

VQE achieved chemical accuracy (within 1 kcal/mol of CCSD(T)) where DFT did not. Three key bottlenecks limit scaling:

• Shot noise: 200-term Pauli Hamiltonians required ~50,000 shots per energy evaluation on IBM hardware
• Circuit depth: UCCSD depth grows steeply with active space size; 10 spatial orbitals already produce thousands of gates
• Barren plateaus: Gradient magnitudes shrink exponentially for larger circuits without specialized initialization

A full kinase binding site involves 50-100 active orbitals, requiring 100-200 logical qubits. Sanofi’s near-term path uses VQE fragment results as training data for ML potentials, improving force field accuracy without requiring full fault-tolerant quantum treatment.

Learn more: Qiskit Reference