⚛️ Quantum Approximate Optimization for Protein Folding

The HP Model — Reducing Protein Folding to a Combinatorial Puzzle

Real protein folding is governed by a bewildering combination of hydrogen bonding, van der Waals forces, electrostatics, backbone torsion strain, and solvent entropy across thousands of atoms — a problem AlphaFold solves with deep learning trained on hundreds of thousands of experimental structures. The HP (hydrophobic-polar) lattice model, introduced by Ken Dill in 1985, strips this down to the single dominant physical driver believed to govern the earliest stages of folding: the hydrophobic effect, in which nonpolar residues cluster together to minimize contact with water.

• Dill, 1985: Model origin (Biochemistry; lattice HP model)
• 20 → 2: Alphabet reduction (amino acids → H or P)
• 2D square / 3D cubic: Lattice type (simplest common choices)
• proven: HP folding is NP-hard (Crescenzi et al. 1998)

Coarse-graining amino acids and defining the HP energy function

The HP model reduction:

• Each of the 20 natural amino acids is classified as either Hydrophobic (H) — e.g., Leu, Ile, Val, Phe, Met — or Polar (P) — e.g., Ser, Thr, Asn, Gln, Lys, Glu, Arg • A protein sequence becomes a binary string, e.g., HPHPPHHPHH for a 10-residue chain • The backbone is placed on a lattice (2D square lattice for pedagogical/small demonstrations; 3D cubic or face-centered-cubic for more realistic studies) as a self-avoiding walk: consecutive residues occupy adjacent lattice sites, and no two residues may occupy the same site

HP energy function:

E = -Σ_{(i,j) non-adjacent in sequence, adjacent on lattice} ε · [h_i=H][h_j=H]

• Every pair of H residues that are NOT sequence-neighbors (i.e., not directly bonded) but ARE spatially adjacent on the lattice (a "topological contact") contributes a favorable energy ε (conventionally ε=-1 per contact) • P-P and H-P contacts contribute zero energy in the simplest version of the model • The folding problem becomes: find the self-avoiding walk conformation that MAXIMIZES the number of non-sequential H-H contacts (equivalently, minimizes E) • This directly encodes the core hydrophobic-collapse intuition: hydrophobic residues want to bury themselves away from solvent by clustering in a compact core, exactly as observed in real globular protein cores

Computational complexity: • Even in this radically simplified form, finding the global energy minimum conformation of the HP model on a 2D or 3D lattice has been formally proven NP-hard (Crescenzi, Goldman, Papadimitriou, Piccolboni, Yannakakis, 1998) — the number of possible self-avoiding walks grows exponentially with chain length (roughly 2.6^N in 2D), so brute-force enumeration is intractable beyond roughly 20-30 residues even on classical supercomputers • This NP-hardness, despite the model's extreme simplicity, is precisely what makes HP folding an attractive target for quantum optimization algorithms like QAOA: it is a genuine combinatorial optimization problem, not merely a toy, even though it captures only a sliver of real protein physics • Crucially, solving the HP model optimally tells you almost nothing directly transferable to real drug-relevant protein structures — it is a testbed for the optimization algorithm, not a competitor to structure prediction tools

Encoding Backbone Conformations as Qubits — Turn-Based and Diamond Encodings

To let a quantum computer search over lattice conformations, every possible fold of the chain must be represented as a bitstring on qubits. The most common approach encodes not residue positions directly, but the sequence of "turns" the backbone takes at each step — a compact representation that scales linearly, rather than exponentially, in the number of qubits needed per residue.

• 2: Qubits per turn (2D) (4 possible directions)
• 2–3: Qubits per turn (3D) (6 possible directions)
• ~2(N−2): Total qubits, N residues (2D square lattice, turn encoding)
• ~8: 6-residue chain qubit count (typical small demo size)

Turn-direction encoding and the diamond encoding for reduced qubit count

Absolute-position encoding (naive, qubit-expensive): • Represent each residue's (x,y) lattice coordinate directly in binary • Requires O(N log L) qubits where L is lattice size — and needs many additional penalty terms to enforce self-avoidance (no two residues at the same coordinate) and chain connectivity (consecutive residues must be lattice-adjacent) • Rapidly becomes qubit-expensive and constraint-heavy; rarely used in practice for quantum approaches

Turn-direction encoding (standard in QAOA protein folding literature, e.g., Robert et al. 2021 "Resource-efficient quantum algorithm for protein folding", npj Quantum Information):

• Rather than encoding absolute positions, encode the direction of each bond relative to the previous one: at each step, the chain can go straight, turn left, or turn right (in a reduced "diamond" lattice representation that eliminates some redundant symmetric conformations) • Each turn requires 2 qubits (4 possible direction states: e.g., 00=straight-continue, 01=left, 10=right, 11=reserved/unused in some encodings) using one-hot or binary encoding schemes • For a chain of N residues, there are N-1 bonds and therefore N-2 independent turn choices (the first bond direction is fixed by convention to remove overall rotational symmetry) — giving roughly 2(N-2) qubits • Example: a 6-residue HP chain needs about 2×4=8 qubits for turn encoding — small enough to simulate exactly or run on today's small NISQ devices • The diamond encoding (used in the Robert et al. paper, and demonstrated on IBM hardware) further removes qubits associated with turns that would immediately violate the lattice geometry, achieving the most qubit-efficient known encoding for this problem class

Constraint terms: • Overlap/self-avoidance penalty: an additional quadratic penalty term added to the cost Hamiltonian penalizes bitstrings whose decoded conformation would place two residues on the same lattice site • These penalty terms are large-coefficient additions to the QUBO cost function (Stage 3) — chosen large enough that no optimal (H-H-contact-maximizing) solution can also violate self-avoidance, but not so large that they dominate the energy landscape and defeat the optimizer • Balancing penalty weight against the true HP energy term is itself a nontrivial hyperparameter tuning problem specific to each problem instance

From HP Energy to a QUBO Cost Function and Ising Hamiltonian

QAOA requires its target problem to be expressed as an Ising Hamiltonian — a sum of Pauli-Z terms (and products of Z terms) acting on qubits, mirroring a classical spin-glass energy function. The HP folding energy plus its self-avoidance penalty terms must therefore be algebraically rewritten first as a Quadratic Unconstrained Binary Optimization (QUBO) problem over 0/1 variables, then converted to the equivalent ±1 spin Ising form via a standard linear substitution.

• {0,1}: QUBO variable domain (binary decision variables)
• {−1,+1}: Ising variable domain (spin variables)
• sᵢ = 1 − 2xᵢ: Substitution (binary → spin mapping)
• H_C = ΣJᵢⱼZᵢZⱼ + Σhᵢ Zᵢ: Cost Hamiltonian form (diagonal in computational basis)

QUBO construction, penalty methods, and the Ising cost Hamiltonian

QUBO formulation of HP folding:

Minimize: E(x) = E_HP(x) + λ · E_penalty(x)

where x is the vector of binary turn-encoding qubits, and:

• E_HP(x): the negative H-H contact count, expressed as a quadratic function of the binary variables — a contact between residue i and residue j exists only if the decoded turn sequence places them adjacent, which itself is a product of several binary turn-indicator variables, expanded into pairwise (quadratic) and sometimes higher-order terms that must themselves be reduced to quadratic form via auxiliary variables

• E_penalty(x): quadratic penalty terms enforcing (a) self-avoidance — no two residues share a lattice site, and (b) one-hot validity — each turn's 2-qubit encoding must represent exactly one of the allowed physical turn directions, not an invalid combination

• λ: penalty weight coefficient, tuned empirically large enough to forbid invalid conformations from ever having lower cost than any valid one

Ising Hamiltonian conversion (required for QAOA):

Substitute each binary variable x_i ∈ {0,1} with a spin variable via x_i = (1 − s_i)/2, where s_i ∈ {−1,+1}. This transforms the QUBO polynomial into:

H_C = Σ_{i<j} J_ij Z_i Z_j + Σ_i h_i Z_i + constant

• Z_i is the Pauli-Z operator on qubit i (eigenvalues +1/-1 correspond to spin up/down, i.e., bit value 0/1) • J_ij (coupling strengths) and h_i (local fields) are determined directly by the QUBO coefficients from the HP energy and penalty terms • This H_C is diagonal in the computational (Z) basis — critically, this means it commutes with itself at all times, and its exponential e^{-iγH_C} can be implemented as a product of single-qubit RZ gates and two-qubit ZZ-interaction gates (implemented via CNOT-RZ-CNOT sandwiches), regardless of circuit depth

This QUBO/Ising formulation is identical in spirit to how QAOA is applied to other combinatorial problems — Max-Cut, portfolio optimization, vehicle routing — the HP folding problem is simply a domain-specific instance of the same general mathematical machinery, which is precisely why QAOA research groups adopted it as an accessible, interpretable "hello world" for quantum protein science.

The QAOA Circuit — Alternating Cost and Mixer Unitaries

QAOA (Quantum Approximate Optimization Algorithm, Farhi, Goldstone, Gutmann 2014) is itself a hybrid variational algorithm structurally similar to VQE, but purpose-built for combinatorial optimization on Ising-form cost Hamiltonians. Starting from an equal superposition over all bitstrings, it applies p alternating rounds of two unitaries — one that imprints phase according to the cost Hamiltonian, one that spreads (mixes) amplitude across bitstrings — engineered so that measurement probability concentrates on low-cost (low-energy, good-fold) bitstrings as p and the tuned parameters improve.

• Farhi et al. 2014: Original QAOA paper (arXiv:1411.4028)
• 2p: Circuit parameters (γ₁..γₚ, β₁..βₚ)
• e^{−iγH_C}: Cost unitary (diagonal phase rotation)
• e^{−iβH_M}, H_M=ΣXᵢ: Mixer unitary (transverse-field mixing)

Circuit structure, the adiabatic connection, and depth-p scaling

QAOA circuit construction:

1. Initial state: |+⟩^⊗n = H^⊗n |0⟩^⊗n — an equal superposition over all 2^n possible turn-encoding bitstrings, i.e., every self-avoiding-walk-attempt (valid or invalid) simultaneously

2. Alternating layers (repeated p times):

Cost unitary: U_C(γ) = e^{-iγH_C} • Implements the Ising cost Hamiltonian as a phase: each computational basis state |x⟩ (each candidate fold) picks up a phase e^{-iγE(x)} proportional to its HP+penalty energy E(x) • Implemented as a circuit of RZ gates (for the h_i local field terms) and CNOT-RZ-CNOT gadgets (for the J_ij ZZ coupling terms) — directly derived from the Ising Hamiltonian structure in Stage 3

Mixer unitary: U_M(β) = e^{-iβH_M}, typically H_M = Σ_i X_i • Implemented as a layer of RX(2β) rotations on every qubit • Its role is to move amplitude between different bitstrings, allowing the algorithm to explore the solution space rather than getting stuck imprinting phase on a fixed superposition

3. Final circuit: |γ,β⟩ = U_M(β_p)U_C(γ_p)...U_M(β_1)U_C(γ_1)|+⟩^⊗n

Connection to the adiabatic theorem: • QAOA is formally a Trotterized approximation of adiabatic quantum computation: as p→∞ with appropriately scheduled (γ,β), QAOA provably converges to the exact ground state of H_C (the global optimum fold), mirroring how a slow adiabatic sweep from a simple Hamiltonian to H_C would remain in the ground state throughout • In practice, p is kept small (1-10) on NISQ hardware because each layer adds substantial two-qubit gate depth, and gate errors compound multiplicatively with depth

Depth-p scaling and protein-folding-specific circuit cost: • Each layer requires implementing every J_ij term in the Ising Hamiltonian as a two-qubit gadget — for an 8-qubit HP folding instance, this is on the order of 10-20 CNOT gates per layer just for the cost unitary, plus the constraint-penalty terms which typically dominate the term count • A p=3 circuit for a 6-residue HP chain can already require 60-100+ two-qubit gates — pushing against the coherence limits of current superconducting and trapped-ion hardware for anything beyond toy chain lengths • This qubit and gate-depth cost is why essentially all published QAOA protein folding demonstrations to date use chains of 4-8 residues, far short of any pharmaceutically informative peptide length

Closing the Loop — Parameter Optimization, Sampling, and the Reality Check Against AlphaFold

Like VQE, QAOA is only half quantum: the parameters (γ,β) controlling the cost and mixer unitaries are tuned by a classical outer-loop optimizer that repeatedly runs the quantum circuit, measures the expected cost, and adjusts parameters to drive that expectation lower — after which many samples are drawn from the tuned circuit and the lowest-energy valid conformation observed becomes the predicted fold.

• COBYLA, Nelder-Mead, Bayesian: Outer-loop optimizers (same family as VQE)
• 4–7 residues: Published demo scale (IBM, D-Wave hybrid studies)
• thousands of residues: AlphaFold2 scale (full proteins, near-experimental accuracy)
• 2021: AlphaFold2 publication (Jumper et al., Nature)

Parameter tuning, sampling readout, and scale comparison with classical structure prediction

Classical optimization of QAOA parameters:

• The expected cost ⟨γ,β|H_C|γ,β⟩ is estimated by repeated circuit execution and measurement, exactly as in VQE's expectation-value estimation step • A classical optimizer (COBYLA is common for the smoothness of small-p QAOA landscapes; Bayesian optimization is increasingly favored because QAOA parameter landscapes can have many local minima) proposes new (γ,β) to reduce ⟨H_C⟩ • Known good initialization strategies exist: parameters tend to vary smoothly and predictably with p, so parameters optimized at depth p can often be interpolated to warm-start optimization at depth p+1 ("INTERP" strategy, Zhou et al. 2020) — substantially reducing the classical optimization burden

Final readout: • Once (γ,β) converge, the tuned circuit is sampled thousands of times (each sample collapses the superposition to one bitstring — one candidate fold) • Each sampled bitstring is decoded back into a lattice conformation; invalid conformations (self-avoidance or turn-encoding violations) are discarded • Among valid conformations, the one with the lowest HP energy (most H-H contacts) observed across all samples is reported as the QAOA-predicted fold • For small chains (≤8 residues), this can be directly checked against the true global optimum found by classical brute-force enumeration — published demonstrations (e.g., Robert et al. 2021 on IBM hardware, and D-Wave quantum annealing HP-folding studies) report successfully recovering the known optimal or near-optimal fold for chains up to about 7 residues

Scale comparison — the honest gap to real structural biology: • A 6-residue HP model instance with 8 qubits is scientifically a proof-of-concept for the optimization method, not for protein structure prediction • AlphaFold2 (Jumper et al., Nature 2021) and AlphaFold3 predict full 3D atomic structures for proteins of hundreds to thousands of residues, at accuracy often competitive with experimental X-ray crystallography (median backbone RMSD often <1-2 Å on CASP benchmark targets), trained on the full Protein Data Bank plus massive multiple-sequence-alignment evolutionary data — a fundamentally different, data-driven paradigm rather than physics-based combinatorial search • Even optimistic extrapolations of QAOA/quantum-annealing hardware roadmaps (hundreds to low thousands of logical, error-corrected qubits) fall dramatically short of the qubit counts a full all-atom, off-lattice conformational search of a 100+ residue protein would require — likely millions of qubits with error correction overhead included, using the most straightforward encodings • The credible near/medium-term role for QAOA-style quantum optimization in structural biology is not competing with AlphaFold on full structure prediction, but as a specialized subroutine for narrow combinatorial sub-problems — e.g., side-chain rotamer optimization, docking pose search, or peptide macrocycle conformer search — where problem sizes are naturally small and classical methods already struggle combinatorially

A useful way to frame QAOA protein folding for a general audience: it demonstrates real, working quantum optimization on a real (if drastically simplified) combinatorial formulation of folding, successfully recovering known-optimal answers for toy problems — a valid and important algorithmic proof-of-concept. It is not, and is not intended by its own authors to be, a competitor to AlphaFold or classical molecular dynamics at any biologically relevant protein scale today.