Roche: Quantum Computing for Protein Folding and Peptide Energy Calculation

The Problem

Proteins fold from linear chains of amino acids into three-dimensional structures, and that structure determines function. Misfolded proteins are implicated in Alzheimer’s disease, Parkinson’s disease, and a range of other conditions. Understanding the folding energy landscape well enough to predict structure from sequence, or to identify misfolding pathways, is one of computational biology’s central challenges.

Classical approaches (AlphaFold, molecular dynamics with force fields like Amber or CHARMM) have made remarkable progress. But they operate on parameterized classical energy functions, not first-principles quantum mechanics. For small peptides where quantum effects in bond angles and torsion potentials matter, more fundamental methods are needed.

Roche’s computational biology team explored two quantum computing approaches at different levels of the problem:

• Coarse-grained folding: D-Wave quantum annealing for lattice protein models
• Fine-grained electronic structure: VQE on IBM Quantum for peptide bond energy calculations

D-Wave for Lattice Protein Folding

The HP (hydrophobic-polar) lattice model is a simplified protein representation. Each residue is classified as H (hydrophobic) or P (polar). The protein folds on a 2D or 3D lattice, and the energy is lowered by H-H contacts between non-adjacent residues. Finding the lowest-energy conformation is NP-hard.

This maps naturally to a QUBO, which D-Wave’s quantum annealer can tackle directly.

import dimod
from dwave.system import DWaveSampler, EmbeddingComposite
import numpy as np

# HP lattice folding for a 6-residue chain: H-P-H-H-P-H
# Encoded as turn directions on a 2D square lattice
# Variables: t_i in {0,1,2,3} = {right, up, left, down} at each step
# Binary encoding: each turn needs 2 bits -> 10 binary variables for 5 turns

sequence = "HPHHPH"
n_residues = len(sequence)
n_turns = n_residues - 1

# Simplified QUBO encoding for 2D square lattice turns
# Turn variables: each turn encoded as 2 bits (t0, t1)
# 00=right, 01=up, 10=left, 11=down

# Build adjacency penalty matrix to reward H-H contacts
# and constraint terms to avoid self-intersection
def build_hp_qubo(sequence):
    n = len(sequence)
    n_turn_bits = 2 * (n - 1)  # 2 bits per turn direction
    Q = {}

    # Constraint: valid lattice path (no self-intersections)
    # Encoded as penalties on incompatible turn combinations
    for i in range(0, n_turn_bits - 2, 2):
        # Penalize U-turns (right then left, or up then down)
        Q[(i, i+2)]   = Q.get((i, i+2), 0) + 2.0
        Q[(i+1, i+3)] = Q.get((i+1, i+3), 0) + 2.0

    # Reward H-H contacts (simplified: reward turns that bring
    # hydrophobic residues adjacent on the lattice)
    for i in range(n):
        for j in range(i + 2, n):
            if sequence[i] == 'H' and sequence[j] == 'H':
                # Approximate: reward when separation is even (same axis)
                if (j - i) % 2 == 0:
                    bit_i = 2 * max(0, i - 1)
                    if bit_i < n_turn_bits:
                        Q[(bit_i, bit_i)] = Q.get((bit_i, bit_i), 0) - 1.0

    return Q

Q = build_hp_qubo(sequence)

# Convert to BinaryQuadraticModel for D-Wave
bqm = dimod.BinaryQuadraticModel.from_qubo(Q)

# On real D-Wave hardware (requires API key and access):
# sampler = EmbeddingComposite(DWaveSampler())
# sampleset = sampler.sample(bqm, num_reads=1000, annealing_time=20)

# Simulate with classical exact solver for validation
from dimod import ExactSolver
sampleset = ExactSolver().sample(bqm)
best = sampleset.first
print(f"Best energy: {best.energy:.4f}")
print(f"Best conformation bits: {list(best.sample.values())}")

# Reconstruct 2D coordinates from turn bits
def decode_path(sample, n_residues):
    directions = [(1,0), (0,1), (-1,0), (0,-1)]  # right, up, left, down
    positions = [(0, 0)]
    for t in range(n_residues - 1):
        b0 = sample.get(2*t, 0)
        b1 = sample.get(2*t + 1, 0)
        d = directions[2*b0 + b1]
        last = positions[-1]
        positions.append((last[0] + d[0], last[1] + d[1]))
    return positions

path = decode_path(best.sample, n_residues)
print(f"Folded path coordinates: {path}")

VQE for Peptide Bond Electronic Structure

At a finer scale, Roche used VQE to compute the electronic energy of short peptide chains as a function of backbone dihedral angles (phi and psi in Ramachandran notation). For a 3-residue chain (alanine tripeptide as a model), the active space involves the peptide bond pi system.

from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock
from qiskit_algorithms.minimum_eigensolvers import VQE
from qiskit_algorithms.optimizers import COBYLA
from qiskit_aer.primitives import Estimator

# Model the peptide bond in formamide (H2N-CHO) as a minimal peptide proxy
# Formamide captures the resonance structure of the peptide bond C-N
driver = PySCFDriver(
    atom="C 0.0 0.0 0.0; O 1.22 0.0 0.0; N -0.68 1.14 0.0; "
         "H -0.3 -1.0 0.0; H -1.68 1.14 0.0; H -0.18 2.05 0.0",
    basis="sto-3g",
    charge=0,
    spin=0,
)
problem = driver.run()

# Active space: 6 electrons in 4 orbitals (pi system + lone pairs)
transformer = ActiveSpaceTransformer(num_electrons=6, num_spatial_orbitals=4)
active_problem = transformer.transform(problem)

mapper = JordanWignerMapper()
hamiltonian = mapper.map(active_problem.second_q_ops()[0])

print(f"Qubit count for peptide bond active space: {hamiltonian.num_qubits}")

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

estimator = Estimator()
optimizer = COBYLA(maxiter=400)
vqe = VQE(estimator=estimator, ansatz=ansatz, optimizer=optimizer)
result = vqe.compute_minimum_eigenvalue(hamiltonian)

vqe_energy = result.eigenvalue.real
print(f"VQE peptide bond energy: {vqe_energy:.6f} Ha")
# DFT/B3LYP reference for formamide pi system: compare in kcal/mol

Comparison with Classical References

Roche validated both approaches against established classical methods:

D-Wave lattice folding vs. exhaustive search: For HP sequences up to 8 residues (where exhaustive search is tractable), D-Wave found the known ground-state conformation in all 12 tested instances using 1000 annealing reads. At 12 residues, D-Wave found the best known solution 9 out of 12 times, matching simulated annealing performance.

VQE vs. DFT (B3LYP/6-31G*): For 3-residue peptide chain fragments, VQE energies agreed with DFT reference values within 0.5 kcal/mol. DFT was run using Gaussian 16. The 0.5 kcal/mol agreement is within the DFT method’s own uncertainty for these systems, meaning VQE results are consistent with the best available classical reference.

Scaling Gap to Biological Relevance

Neither quantum approach is yet competitive with classical methods at biologically meaningful scale:

• A biologically relevant protein has 100-1000 residues. The D-Wave lattice model at that scale requires thousands of variables, demanding embedding that exceeds current D-Wave chip connectivity.
• Fine-grained VQE for a full 3-residue alanine tripeptide (not just the pi system proxy) requires 40-60 qubits in the active space. Current hardware noise budgets limit useful VQE to around 10-15 qubits.

Roche treats current results as protocol validation. The quantum workflows, benchmarking infrastructure, and QUBO formulations developed now will be directly reusable when hardware improves.

Results

Key findings from Roche’s dual-platform quantum program:

• D-Wave lattice folding: matched exhaustive search on all HP instances up to 8 residues using 1000 annealing reads
• VQE peptide bond calculations: 0.5 kcal/mol agreement with DFT for 3-residue chains
• Computational biology team published benchmark protocols usable by future research groups
• Both approaches remain far from biological scale; quantum hardware needs 2-3 generations of improvement
• Roche positions quantum computing as a long-range investment in computational biology infrastructure

Two-Platform Strategy

Roche’s use of both gate-based (IBM Quantum) and annealing (D-Wave) hardware reflects the different problem structures:

• D-Wave (annealing): suited for combinatorial search over discrete conformational space. The HP lattice model is a natural QUBO. Annealing does not require deep circuits and is less sensitive to decoherence than gate-based methods.
• IBM Quantum (gate-based): suited for electronic structure calculations where the continuous energy landscape requires VQE-style variational optimization. Provides a path to fault-tolerant quantum chemistry as hardware improves.

This bifurcated strategy lets Roche extract maximum value from each hardware paradigm while tracking progress across both.

Learn more: Qiskit Reference | D-Wave Ocean SDK Reference | VQE Algorithm Guide