Resource-Efficient Quantum Algorithm for Protein Folding
Anton Robert, Panagiotis Kl. Barkoutsos, Stefan Woerner, Ivano Tavernelli
RResource-Efficient Quantum Algorithm for Protein Folding
Anton Robert,
1, 2
Panagiotis Kl. Barkoutsos, Stefan Woerner, and Ivano Tavernelli ∗ IBM Research GmbH, Zurich Research Laboratory, Säumerstrasse 4, 8803 Rüschlikon, Switzerland PASTEUR, Département de chimie, École Normale Supérieure,PSL University, Sorbonne Université, CNRS, 75005 Paris, France (Dated: August 7, 2019)
Predicting the three-dimensional (3D) structureof a protein from its primary sequence of aminoacids is known as the protein folding (PF) prob-lem. Due to the central role of proteins’ 3D struc-tures in chemistry, biology and medicine applica-tions (e.g., in drug discovery) this subject has beenintensively studied for over half a century.
Al-though classical algorithms provide practical so-lutions, sampling the conformation space of smallproteins, they cannot tackle the intrinsic NP-hardcomplexity of the problem, even reduced to itssimplest Hydrophobic-Polar model. While fault-tolerant quantum computers are still beyond reachfor state-of-the-art quantum technologies, there isevidence that quantum algorithms can be success-fully used on Noisy Intermediate-Scale Quantum(NISQ) computers to accelerate energy opti-mization in frustrated systems.
In this work,we present a model Hamiltonian with O ( N ) scal-ing and a corresponding quantum variational algo-rithm for the folding of a polymer chain with N monomers on a tetrahedral lattice. The model re-flects many physico-chemical properties of the pro-tein, reducing the gap between coarse-grained rep-resentations and mere lattice models. We use a ro-bust and versatile optimisation scheme, bringing to-gether variational quantum algorithms specificallyadapted to classical cost functions and evolution-ary strategies (genetic algorithms), to simulate thefolding of the 10 amino acid Angiotensin peptideon 22 qubits. The same method is also successfullyapplied to the study of the folding of a 7 amino acidneuropeptide using 9 qubits on an IBM Q 20-qubitquantum computer. Bringing together recent ad-vances in building gate-based quantum computerswith noise-tolerant hybrid quantum-classical algo-rithms, this work paves the way towards accessibleand relevant scientific experiments on real quantumprocessors. The solution of Levinthal’s paradox through a biassearch in configuration space demands a very fine de-scription of the interactions in the cellular environmentto correctly drive the search in the rugged energy land-scape of a protein. GPU-assisted sampling methodsof well parametrised coarse-grained models can provide use-ful insights regarding the protein’s native conformation,but folding a small protein in state-of-the-art simulationscomes at a very high computational cost.
Protein latticemodels reduce the conformational space and obviate thehigh computational cost of an off-lattice exhaustive sam-pling.
Quantum algorithms cannot ignore those simpli-fications because of curently available quantum resources. Perdomo-Ortiz et al. paved the way towards the construc-tion of spin Hamiltonian to find the on-lattice heteropoly-mer’s low-energy conformations using quantum devices, butwith unattainable high costs for NISQ computers.
Re-cently, the amount of resources needed for the simulationof polymer lattice models was reduced, however still main-taining an exponential cost in terms of number of qubitsand gates needed.
These methods were used to fold acoarse-grained protein model with and amino acid se-quences on a 2D and 3D lattice, respectively, using a quan-tum annealer. These experiments required and qubits and led to a final population of . and . for the corresponding ground state structures, using di-vide and conquer strategies. More recently, Fingerhuth andcoworkers proposed another approach based on the Quan-tum Approximate Optimization Algorithm (QAOA) us-ing a problem-specific alternating operator ansatz to modelprotein folding. Employing the same model proposed in they succeeded in folding a amino acid protein model ona 2D square lattice.In this work, we present a coarse-grained model forprotein folding which is suited to the representation ofbranched heteropolymers comprised of N monomers on atetrahedral (or “diamond”) lattice. This choice is moti-vated by the chemical plausibility of the angles enforcedby the lattice ( . o for bond angles, o or o for di-hedrals), which allows an all-atom description for a widerange of chemical and biological compounds. Given themodest resources of actual quantum devices, a two-centeredcoarse-grained description of amino acids (backbone andside chain) was used to mimic the protein sequence. Everymonomer is depicted by one or multiple beads that can havea defined number of ‘color shades’ corresponding to differ-ent physical properties like hydrophobicity and charge. The configuration qubits.
As for the previous models inliterature, a polymer configuration is grown on the latticeby adding the different beads one after the other and encod-ing, in the qubit register the different “turn” t i that definesthe position of the bead i + 1 relatively to the previousbead i . Using a tetrahedral lattice, we distinguish two setsof nonequivalent lattice points A and B (see Fig. 1). Atthe A sites, the polymer can only grow along the directions t i ∈ { , , , } while at site B the possible directions are t i ∈ { ¯0 , ¯1 , ¯2 , ¯3 } . Along the sequence, the A and B sitesare alternated so that we can use the convention that A (respectively B ) sites correspond to even (odd) i . Withoutloss of generality, the first two turns can be set to t = ¯1 and t = 0 due to symmetry degeneracy. To encode theturns, we assign one qubit per axis t i = q i − q i − q i − q i (Fig. 1(c)). Therefore, the total number of qubits requiredto encode a conformation q cf corresponds to N cf = 4( N − .If the monomers are described by more than one bead, the a r X i v : . [ qu a n t - ph ] A ug N N u m b e r o f P a u li s t r i ng s N t a ( N b ) a ( N b ) N N u m b e r o f qub i t s sparsedense A
Tetrahedral lattice polymer model. (a)
La-belling of the coordinate systems at the sub-lattices A and B . (b) Typical polymer conformation (10 monomers). The red anddark green dashed lines represent a subset of inter-bead inter-actions considered in our model. Side chain beads are shown inorange. (c)
Example of turn encoding. (d)
Number of qubitsrequired by the sparse (3-local terms) model (blue) and its dense(5-local terms) variant (orange) as a function of the number ofmonomers. (e)
Number of Pauli strings with respect to thenumber of monomers for the sparse and dense encoding mod-els. The parameters of the fit are ( a, b ) = (0 . , . . Theexponential curve (black dotted) is given as a reference. same formula holds by replacing N with the total numberof beads in the polymer. A denser encoding of the polymerchain using only N − configuration qubits is presentedin the SI. The interaction qubits.
To describe the interactions, weintroduce a new qubit register q in , composed of q ( l ) i,j for each l th nearest neighbour ( l -NN) interaction on the lattice (seered and green dashed lines for l = 1 and l = 2 in Fig. 1(b))between beads i and j . The use of these registers will beexplained in connection to the definition of the interactionenergy terms. The number of qubits constituting the inter-action register, N in , is entirely determined by the skeletonof the polymer (i.e. including the side chains), regardlessof the beads’ color, and scales as O ( N ) . Note that two -NN beads occupy positions on different sub-lattices ( A or B ). On the other hand, for l > all beads of both sub-lattices can potentially interact. Given a primary sequence,the pairwise interaction energies (cid:15) ( l ) i,j between the beads atdistance l can be arbitrarily defined to reproduce a fold ofinterest or it can be adapted from pre-existing models, likethe one proposed by Miyazawa and Jernigan (MJ) for -NNinteractions. The Hamiltonian.
The next step defines the qubitHamiltonian that describes the energy of a given fold de-fined by the sequence of beads (fixed) and the encodedturns. Penalty terms are applied when physical constraintsare violated (e.g., when beads occupy the same position onthe lattice), and physical interactions (attractive or repul-sive in nature) are applied when two beads occupy neigh-bouring sites or are at distance l > , where l is the num-ber of NN in real space. The different contributions to thepolymer Hamiltonian are therefore (with q = { q cf , q in }), H ( q ) = H gc ( q cf ) + H ch ( q cf ) + H in ( q ) . (1)The definitions of the geometrical constraint ( H gc , whichgoverns the growth of the primary sequence with no bifur-cation) and the chirality constraint ( H ch , which enforcesthe correct stereochemistry of the side-chains if present)are given in SI. The interaction energy terms.
For each bead i alongthe sequence the distance to the other beads j (cid:54) = i canuniquely be determined by the state of the N cf configura-tion qubits. To this end, for each pair of beads ( i, j ) we in-troduce a four-dimensional vector (Eq. SI-13), the norm ofwhich uniquely encodes they reciprocal distance d ( i, j ) . Asan example, we consider the energy contributions for -NNinteractions. For each pair of beads ( i, j ) an energy contri-bution of (cid:15) ( l ) ij is added to H i,j in when the distance d ( i, j ) = l .However, a contribution of the form (cid:15) ( l ) ij δ ( d ( i, j ) − l ) can-not be efficiently implemented as a qubit string Hamilto-nian (here δ ( . ) stands for the Dirac delta function). Usingthe set of contact qubits q ( l ) i,j we therefore define an energyterm of the form q ( l ) i,j ( (cid:15) ( l ) ij + λ ( d ( i, j ) − l )) for each value of l and λ (cid:29) (cid:15) ( l ) ij . This definition implies that the contribution (cid:15) ( l ) ij for the formation of the “interaction” ( i, j ) at distance l is only assigned when the contact qubit q ( l ) i,j = 1 and d ( i, j ) = l , simultaneously. For q ( l ) i,j = 1 and d ( i, j ) (cid:54) = l thefactor λ adds a large positive energy contribution that over-comes the stabilizing energy (cid:15) ( l ) ij (the case of d ( i, j ) < l isdetailed in the SI). Finally, in our model we prevent the si-multaneous occupation of a single lattice site by two beads,as discussed in the SI. In a nutshell, we only prevent over-laps that occur in the vicinity of an interaction pair. If q ( l ) i,j = 1 , we apply penalty functions so that i and j + 1 cannot overlap when l = 1 , for instance. The folding algorithm.
The solution to the folding prob-lem is the ground state of the Hamiltonian H ( q ) and there-fore lies in the N cf dimensional space of the configurationqubits. To find this solution we prepare a variational cir-cuit, comprising both the configurational and the interac-tion registers, which is composed by an initialization blockwith Hadamard gates and parametrized single qubit R Y gates followed by an entangling block and another set ofsingle qubit rotations. We denote by θ = ( θ cf , θ in ) theset of angles of size n where n = N cf + N ct is the totalnumber of qubits. Differently to the quantum mechani-cal case, for the solution of the ‘classical problem’ (e.g.,folding) we do not need an estimate of the Hamiltonian ex-pectation value, but we only require the sampling of thelow energy tail of the energy distribution. Therefore, the -9.88 -9.69 -8.18 -6.48 -5.84 -5.1 -3.52 -2.32 -1.7Energy of conformation020406080 P o pu l a t i o n % H P F H L D -2.32 -3.48 -3.4R -1.7 -2.16V -6.29 -6.48Y -3.52V -6.48 (a) (c)(b) ECVaR ( ! ) createnextgenerationreplace parent if offspringsbetter | i
QuantumClassical p
Starting from a random population (up-center) of circuitparameters { θ } , every parent, θ p , undergoes a parametrized recombination with other individuals according to the proceduresdetailed in Section Materials and Methods. The corresponding trial wavefunctions are generated by the quantum circuit as describedin the main text and measured to estimate the new CVaR objective function. They constitute the selection criteria of whether toreplace a parent by its offspring for the new generation. (b) Folding of the ten amino acid Angiotensin peptide. Energydistribution at the convergence of the low-energy folds for the population obtained with the CVaR-VQE algorithm and the DEoptimizer. The results were obtained using (blue) and measurements (orange). Simulations were carried out using arealistic parametrization of the noise. The binary strings ( q , q , q , q , q , q , q , q , q , ) associated with the differentbars represent the contact qubits (see main text) that entirely define the conformation energies. The numbers labelling the barscorrespond to the exact degeneracy of the conformations. The total percentage of finding low-energy conformations (energy below0) adds up to . (small sampling, blue) and (large sampling, orange). The fittest individual in the population collapses tothe ground state with a probability of . % (see SI). (c) Primary sequence of Angiotensin. Each amino acid is assigned a differentcolour to match its specific physical properties. The letters stand for Aspartic-Acid (D), Arginine (R), Valine (V), Tyrosine (Y),Histidine (H), Proline (P), Phenylalanine (F), and Leucine (L). (d)
Pairwise interaction matrix for Angiotensin constructed usingthe MJ model (Table 3 in ). optimization of the angles θ is performed using a modifiedversion of the Variational Quantum Eigensolver (VQE) algorithm named Conditional Value-at-Risk (CVaR) VQEor simply CVaR-VQE. Briefly, CVaR defines an objectivefunction based on the average over the tail of a distributiondelimited by a value α (see histogram in Fig. 2(a)) whichis denoted CVaR α ( θ ) = (cid:104) ψ ( θ ) | H ( q ) | ψ ( θ ) (cid:105) α . Compared toconventional VQE, CVaR-VQE provides a drastic speed-up to the optimization of diagonal Hamiltonians as shownin. The classical optimization of the gate parameters isperformed using a Differential Evolution (DE) optimizer, which mimics natural selection in the space of the angles θ . The optimisation procedure is summarized in Fig. 2(a).Note that at each step of the optimization, the wavefunc-tions | ψ ( θ p ) (cid:105) corresponding to the different individuals θ p (Fig. 2(a)) are collapsed during measurement leading tobinary strings, which are uniquely mapped to the corre-sponding configurations and energies. We denote by P f ( p ) the probability for the individual p to find the f th lowestenergy fold at convergence. Scaling.
We define the scaling of the algorithm as thenumber of terms (or Pauli strings), in the n -qubit Hamil-tonian H ( q ) (see also Table I of SI). H ( q ) = N t (cid:88) γ h γ n (cid:79) i =1 q γ i i (2)where h γ real coefficients, q i = (1 − σ zi ) / where σ zi is Z Pauli matrix with γ i ∈ { , } , and N t is the total num-ber of terms. A thorough investigation of the scaling (seeSI) reveals that the geometrical constraints imposed by the tetrahedral lattice give rise to all possible -local termswithin the N cf conformation qubits. Due to the coupling(entanglement) with the interaction qubits the Hamilto-nian locality (i.e. the maximum number of Pauli opera-tors different from the identity in H ( q ) ) is strictly forthe -NN interaction. Moreover, the scaling is bound by N t ∼ N in (cid:0) N cf (cid:1) = O ( N ) even for l -NN interactions, with l ≥ . Fig. 1(d) and (e) respectively report the scaling ofthe proposed model and its qubits requirements. Applications.
We first apply our quantum algorithm tothe simulation of the folding of the 10 amino acid peptideAngiotensin. Using our coarse grained model on the tetra-hedral lattice the simulation of this system would require35 qubits, which is computationally intensive. We thereforeintroduced a denser encoding of the polymer configurationthat requires only 2 qubits per turn t i = q i − q i , reducingthe total number of qubits to 22. This variant generated -local (instead of -local) terms in the qubit Hamiltonianwhile keeping the total number of Pauli strings within anaffordable range for small instances (see Fig. 1(d)). To fur-ther reduce the number of qubits we also integrate the sidechains with the corresponding bead along the primary se-quence and neglect interactions with l > . Each bin ofthe histogram in Fig. 2(b) counts, over the population, theoccurrence of P f ( p ) for the minimum energy fold ( f = 0 )and the next 18 folds (histogram bars) given n s = 128 (blue bars) and n s = 1024 (orange bars) measurementsof the wavefunctions during the minimization. More than of the individuals in the final population can gener-ate the minimal conformations after generations (orangebars), which occurs with a probability max p P f ( p ) = 42 . . C V a R ↵ AverageAll
HHHHHHHHH R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i c . . . . . P r o b a b ili t i e s HHHHHHHHH R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i c HHHHHHHHH R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y R Y q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i q : | i c q q q q q q q q q q (a) (d)(b)(c) . . . . . . . . P r o b a b ili t i e s Best (e)
FIG. 3.
Experiment results for the folding of the 7 amino acid neuropeptide. (a)
Parametrized quantum circuit for thegeneration of the configuration space. The optimal set of qubit gate rotations is used to reconstruct the optimal fold. (b)
Sectionof the 20 qubits IBM Q Poughkeepsie device used in this experiment. Qubit 7 is used to close the loop by swapping with qubit8. (c)
Maximal converged ground state probability for the ground state fold in the population, max p P ( p ) . (d) Evolution of theground state energy during the CVaR-VQE minimization (i.e., number of generations) with α parameter set to . (e) Evolutionof the mean probability (over the population ensemble, (cid:104) P ( p ) (cid:105) ) and of the best individual probability for the ground state fold, max p P ( p ) . The evolution of the percentage throughout the minimiza-tion can be found in the SI. By reducing the number ofmeasurements to 128 shots, we obtained a broader spec-trum of low energy conformations, which still includes theglobal minimum but with a lower probability. Among thelow-energy conformations (with energies below 0), we canclearly identify the formation of an α -helix and a β -sheet(conformations marked with a grey arrow in Fig. 2(b)). Bytuning the interaction matrix (see Fig.S2 in SI), we canfoster the formation of secondary structural elements. The22-qubit Angiotensin system is still too large for encodingin state-of-the-art quantum hardware. To this end we useda smaller 7 amino acid neuropeptide with sequence APRL-RFY (one-letter coding) that can be mapped to 9 qubits.The corresponding CVaR-VQE circuit is shown in Fig. 3(a).As the entangling block we used a closed-loop of CNOTgates that fits the hardware connectivity of the 20-qubitIBM Q device Poughkeepsie (Fig. 3(b)). The mean CVaR α energy value of the population as a function of the numberof generations shows a robust and smooth convergence to-wards the optimal fold (Fig. 3(d)). More importantly, theaverage probability (Fig. 3(e)) of the ground state averagedover the entire population, (cid:104) P ( p ) (cid:105) , increases monotonicallyreaching a final value larger than 20% and with max p P ( p ) peaking up at (see Fig. 3(c) and (e)). Discussion and conclusions.
In this work, we intro-duced a quantum algorithm for the solution of the PFproblem on a regular tetrahedral lattice. The model Hamil-tonain describes a primary coarse-grained protein sequencewhere each beads represents an amino acid. Side chains canalso be modeled by means of an additional bead linked tothe main chain. The interaction between the amino acids(backbone and side chains) can be extended to l -NN (with l > ) along the lattice edges. This enables the model-ing of sophisticated coarse-grained models accounting forLennard-Jones and Coulombic like interactions. We showhow the model can correctly reproduce secondary structureelements through the simple adaptation of the imputed contact map. The number of qubits scales quadraticallywith the number of amino acids , while the number of ele-ments in the Hamltonian scales in O ( N ) . This implies theuse of an unconventional treatment of the overlaps, whichare avoided through the addition of penalty terms. Eventhough the PF problem is a classical optimization prob-lem, the variational quantum algorithm used, CVaR-VQE,drastically reduces the number of measurements required tominimize the classical cost function (instead of the quantummechanical average) and may lead to quantum advantagethrough the use of entanglement. The construction of spe-cific mixing ansatz can drastically speed-up the search inthe configuration space even when the ground state is notentangled but classical. The direct connection betweenqubits and physical properties (configuration and contactsor interactions) allows a rationalisation of initialisation ofthe qubit states and their entanglement, beyond the simplescheme adopted in this preliminary investigation.The locality of the Hamiltonian combined with thefavourable scaling of the qubit resources and the circuitdepth with the number of monomers, make our model thecandidate of choice for the solution of the PF problem onNISQ devices and other quantum technologies. In this workwe performed the simulation of the folding of the 10 aminoacids protein Angiotensin using a realistic model for thenoise of the one and two qubit gate operations. Further-more we used a 20 qubit IBM Q processor to compute thefolding of a 7 amino acid peptide on 9 qubits, which is toour knowledge the largest folding calculation on a NISQdevice using a variational algorithm. The success of thiscalculation demonstrates the potential of our folding algo-rithm and opens up new interesting avenues for the useof quantum computers in the optimization of classical costfunctions using the CVaR-VQE approach combined with agenetic algorithm for the selection of the best fitting varia-tional parameters.
AUTHOR’S CONTRIBUTIONS
A.R., S.W. and I.T. designed the project. A.R. and P.B.performed the experiments and the simulations. All au-thors contributed to the analysis of the results and to thewriting of the manuscript.
ACKNOWLEDGEMENTS
I.T. and P.K.B. acknowledge financial support from theSwiss National Science Foundation (SNF) through thegrant No. 200021-179312. A. R. is much obliged toVladimir Nikolaevitch Smirnov who financially supportedhis work. All authors would like to acknowledge supportfrom the IBM Q network and thank the qiskit developmentteam for discussions regarding the development of the soft-ware. IBM, IBM Q, Qiskit are trademarks of InternationalBusiness Machines Corporation, registered in many juris-dictions worldwide. Other product or service names maybe trademarks or service marks of IBM or other companies.
MATERIALS AND METHODS
The noisy simulations were conducted with α = 1% (resp. α = 0 . ) for the shots (resp. shots) sim-ulation with an “all-to-all“ entangling scheme on Qiskit. All circuits were constructed with a VQE depth of m = 2 .Given a run with n qubits, the size of the population forthe evolutionary algorithm was set to P = 5 mn , a typi-cal size according to the literature. The selection strategyof the DE algorithm is practically identical to the original“current-to-best/1/bin”. ∗ [email protected] M. Levitt and A. Warshel, Nature , 694 (1975). R. Zwanzig, A. Szabo, and B. Bagchi, Proceedings of theNational Academy of Sciences , 20 (1992). D. A. Hinds and M. Levitt, Proceedings of the NationalAcademy of Sciences , 2536 (1992). E. I. Shakhnovich, Physical Review Letters , 2618 (1995). J. N. Onuchic and P. G. Wolynes, Current Opinion in Struc-tural Biology , 70 (2004). K. Lindorff-Larsen, S. Piana, R. O. Dror, and D. E. Shaw,Science , 517 (2011). C. Hyeon and D. Thirumalai, Nature Communications , 487(2011). R. Unger and J. Moult, Bulletin of Mathematical Biology ,1183 (1993). B. Berger and T. Leighton, Journal of Computational Biology , 27 (1998). N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow,A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gam-betta, M. Ganzhorn, A. Kandala, A. Mezzacapo, P. MÃijller,W. Riess, G. Salis, J. Smolin, I. Tavernelli, and K. Temme,Quantum Science and Technology , 030503 (2018). J. Preskill, arXiv:1801.00862 [cond-mat, physics:quant-ph](2018), arXiv: 1801.00862. S. Mandrà and H. G. Katzgraber, Quantum Science andTechnology , 04LT01 (2018), arXiv: 1711.01368. J. King, S. Yarkoni, J. Raymond, I. Ozfidan, A. D.King, M. M. Nevisi, J. P. Hilton, and C. C. McGeoch,arXiv:1701.04579 [quant-ph] (2017), arXiv: 1701.04579. J. King, S. Yarkoni, J. Raymond, I. Ozfidan, A. D.King, M. M. Nevisi, J. P. Hilton, and C. C. Mc-Geoch, arXiv:1701.04579 [quant-ph] (2017), 10.1103/Phys-RevX.6.031015, arXiv: 1701.04579. M. Streif and M. Leib, arXiv:1901.01903 [quant-ph] (2019),arXiv: 1901.01903. C. Levinthal, Journal de Chimie Physique , 44 (1968). J. N. Onuchic, P. G. Wolynes, Z. Luthey-Schulten, and N. D.Socci, Proceedings of the National Academy of Sciences ,3626 (1995). S. Piana, J. L. Klepeis, and D. E. Shaw, Current Opinion inStructural Biology , 98 (2014). Y. Duan, Science , 740 (1998). R. Babbush, A. Perdomo-Ortiz, B. O’Gorman, W. Macready,and A. Aspuru-Guzik, arXiv preprint arXiv:1211.3422 (2012). A. Perdomo Ortiz, C. Truncik, I. Tubert-Brohman, G. Rose,and A. Aspuru-Guzik, Physical Review A (2008). A. Perdomo-Ortiz, N. Dickson, M. Drew-Brook, G. Rose,and A. Aspuru-Guzik, Scientific Reports (2012),10.1038/srep00571. T. Babej, C. Ing, and M. Fingerhuth, arXiv:1811.00713[quant-ph] (2018), 00002 arXiv: 1811.00713. J. Cai, W. G. Macready, and A. Roy, arXiv:1406.2741[quant-ph] (2014), arXiv: 1406.2741. E. Farhi and A. W. Harrow, arXiv:1602.07674 [quant-ph](2016). M. Fingerhuth, T. Babej, and C. Ing, arXiv:1810.13411[quant-ph] (2018), 00002 arXiv: 1810.13411. S. Miyazawa and R. L. Jernigan, Journal of molecular biology , 623 (1996). A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q.Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, NatCommun (2014), 10.1038/ncomms5213. J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. , 023023 (2016). P. K. Barkoutsos, G. Nannicini, A. Robert, I. Tavernelli, andS. Woerner, arXiv:1907.04769 [quant-ph] (2019), 00000 arXiv:1907.04769. R. Storn and K. Price, Journal of Global Optimization ,341 (1997). S. Hadfield, Z. Wang, B. O’Gorman, E. G. Rieffel, D. Ven-turelli, and R. Biswas, Algorithms , 34 (2019), arXiv:1709.03489. G. Aleksandrowicz, T. Alexander, P. Barkoutsos, L. Bello,Y. Ben-Haim, D. Bucher, F. J. Cabrera-Hernández,J. Carballo-Franquis, A. Chen, C.-F. Chen, J. M. Chow,A. D. Córcoles-Gonzales, A. J. Cross, A. Cross, J. Cruz-Benito, C. Culver, S. D. L. P. González, E. D. L. Torre,D. Ding, E. Dumitrescu, I. Duran, P. Eendebak, M. Everitt,I. F. Sertage, A. Frisch, A. Fuhrer, J. Gambetta, B. G.Gago, J. Gomez-Mosquera, D. Greenberg, I. Hamamura,V. Havlicek, J. Hellmers, Ł. Herok, H. Horii, S. Hu,T. Imamichi, T. Itoko, A. Javadi-Abhari, N. Kanazawa,A. Karazeev, K. Krsulich, P. Liu, Y. Luh, Y. Maeng, M. Mar-ques, F. J. Martín-Fernández, D. T. McClure, D. McKay,S. Meesala, A. Mezzacapo, N. Moll, D. M. Rodríguez, G. Nan- nicini, P. Nation, P. Ollitrault, L. J. O’Riordan, H. Paik,J. Pérez, A. Phan, M. Pistoia, V. Prutyanov, M. Reuter,J. Rice, A. R. Davila, R. H. P. Rudy, M. Ryu, N. Sathaye,C. Schnabel, E. Schoute, K. Setia, Y. Shi, A. Silva, Y. Sir-aichi, S. Sivarajah, J. A. Smolin, M. Soeken, H. Takahashi,I. Tavernelli, C. Taylor, P. Taylour, K. Trabing, M. Treinish, W. Turner, D. Vogt-Lee, C. Vuillot, J. A. Wildstrom, J. Wil-son, E. Winston, C. Wood, S. Wood, S. Wörner, I. Y. Akhal-waya, and C. Zoufal, “Qiskit: An open-source framework forquantum computing,” (2019). S. Das, S. S. Mullick, and P. Suganthan, Swarm and Evolu-tionary Computation , 1 (2016). SUPPORTING INFORMATION FOR “RESOURCE-EFFICIENT QUANTUM ALGORITHM FOR PROTEINFOLDING”I. LATTICE MODEL
We call i the backbone bead of index i ∈ { , . . . , N } in the primary sequence of the polymer and denote by i (1) , i (2) , ..., i ( s ) the beads constituting its side chain. Without additional information, discussions about the locality of the Hamiltonianare made with respect to the sparser encoding. A. Conformation encoding
Sparser Encoding
The polymer sequence is generated by specifying the series of turns on the lattice, starting from bead i = 1 (see Fig. 1 of the main text). Each turn t i is encoded on qubits q i − q i − q i − q i , each one representing a directionon the tetrahedral lattice (see Fig. 1 of the main text). One and only one qubit of the four will have value one (while theothers will be set to zero). Accordingly, we encode the turns of the side chain in i as q (0) i (1) q (1) i (1) q (2) i (1) q (3) i (1) , q (0) i (2) q (1) i (2) q (2) i (2) q (3) i (2) ... where the upper index in parenthesis labels the 4 qubits describing the turn of side chain bead i ( k ) and the lower indicesof the form k ( l ) label the bead l along the side chain at the main chain position k . Without loss of generality, we choosethe first two turns to be and . A string of bits defining the entire conformation will have the general form ofEq. SI-1 where each ’residue’ (main and side chain beads) is embraced within square brackets and the side-chain qubitsare within parentheses. Note that the first and the last bead of the main chain (terminal groups) have no side chains.This encoding requires N − − N − qubits to totally define a conformation. [0010] (cid:104) (1)2 (1) q (2)2 (1) q (3)2 (1) q (4)2 (1) ... ) (cid:105) (cid:104) q q q q (q (1)3 (1) q (2)3 (1) q (3)3 (1) q (4)3 (1) ... ) (cid:105) ... (cid:2) q N − − q N − − q N − − q N − (cid:3) (SI-1) Denser Encoding
The turn t i is encoded on two qubits q i − q i . Accordingly, we will encode the turns in a given sidechain by q (1) i (1) q (2) i (1) , q (1) i (2) q (2) i (2) , ... Without loss of generality we choose the first two turns to be and . If the bead onthe main chain does not bear a side chain, another qubit can be saved without breaking any symmetry ( q = 1 ). Theside chains can be encoded with the same convention as for the sparser encoding. A string of bits defining the entireconformation will have the general form of Eq. SI-2 where the side chain qubits are in parentheses. This encoding requires N − − N − qubits to totally define a conformation. [01] (cid:104) (1)2 (1) q (2)2 (1) ... ) (cid:105) (cid:104) q q (q (1)3 (1) q (2)3 (1) ... ) (cid:105) (cid:104) q q (q (1)4 (1) q (2)4 (1) ... ) (cid:105) ... (cid:2) q N − − q N − (cid:3) (SI-2) Turn indicator
It is convenient to introduce an indicator for the axis chosen at turn t i . We will denote f a (q i − , q i ) = f a ( i ) the function that returns if the axis a ∈ { , , , } (see main text, Fig. 1) is chosen at turn t i . For the denserencoding this function is given by f ( i ) = (1 − q i − )(1 − q i ) (SI-3) f ( i ) = q i (q i − q i − ) (SI-4) f ( i ) = q i − (q i − − q i ) (SI-5) f ( i ) = q i − q i . (SI-6)For the sparser encoding, the expression of f a ( i ) is trivial: f ( i ) = q i − (SI-7) f ( i ) = q i − (SI-8) f ( i ) = q i − (SI-9) f ( i ) = q i (SI-10) Distances
The shortest spatial distance between two beads is measured as a function of the number of turns separatingthem along the main chain (see below for the case where the distance is computed between two beads of the side chains).We define as n a ( i, j ) (resp. n ¯ a ( i, j ) ) the number of occurrence of a ∈ { , , , } (reps. ¯ a ∈ { ¯0 , ¯1 , ¯2 , ¯3 } ) along the sequenceseparating i and j with j > i . All intermediate beads will be labeled by k = i, . . . , j − . We assume that the polymerstarts with a bead on the sub-lattice B . Due to the alternation of the two sub-lattices, all sites of the sub-lattice A ( B ) willbe labelled with even (odd) integers. Let’s denote ∆ n a ( i, j ) = n a ( i, j ) − n ¯ a ( i, j ) with n a ( i, j ) (resp. n ¯ a ( i, j ) ) the numberof occurrence of a (reps. ¯ a ) in the sequence between bead i and j . We then have ∆ n a ( i, j ) = j − (cid:88) k = i ( − k f a ( k ) (SI-11)The distance between beads that are placed on side chains, say i ( s ) ( s th bead of side chain at position i of the main chain)and j ( p ) ( p th bead of side chain at position j ) with j > i is obtained by first considering the distance between i and j andthen add or remove the contributions of side chain qubits. Because j > i , the side chain qubits of j (resp. i ) are takeninto account with a plus (resp. minus) sign. Moreover, the direction is defined by the parity of i and j . The generalizedexpression for ∆ n a is therefore ∆ n a ( i ( s ) , j ( p ) ) = ∆ n a ( i, j ) + p (cid:88) l =1 ( − l f a ( j ( l ) ) − s (cid:88) m =1 ( − m f a ( i ( m ) ) (SI-12)For the calculation of the actual distance, is is convenient to introduce the four dimensional vector x ( i, j ) defined as x ( i, j ) = ∆ n ( i, j )∆ n ( i, j )∆ n ( i, j )∆ n ( i, j ) (SI-13)such that d ( i, j ) = || x ( i, j ) || = (cid:80) a ∆ n a ( i, j ) . It is crucial to note that on the tetrahedral lattice there is a bijective mapbetween the values of d ( i, j ) and the Euclidean (through-space) distances r ij in lattice bond units, d ( i, j ) = 0 ⇒ r ij = 0 d ( i, j ) = 1 ⇒ r ij = 1 d ( i, j ) = 2 ⇒ r ij = 2 (cid:113) (cid:39) . d ( i, j ) = 3 ⇒ r ij = (cid:113) (cid:39) . d ( i, j ) = 4 ⇒ r ij = √ (cid:39) . d ( i, j ) = 5 ⇒ r ij = (cid:113) (cid:39) . ... (SI-14) B. Construction of local Hamiltonians
Sparser encoding Hamiltonian
For the sparser encoding, we need to impose that one and only one of the four qubitsthat define a turn is equal to one. On a quantum circuit, this can be achieved easily by using a valid initialization of thequbits together with gate operations that conserved the number of 1’s in each set of 4 qubits that encodes a turn. Whenthis is not possible, a penalty function as show in Eq SI-15 with a large positive λ can be used to impose this constraint. H optional = N − (cid:88) i =3 λ (q i − + q i − + q i − + q i − (SI-15) Growth constraint.
Several constraints are needed to prevent the growth of the chain towards unphysical geometries(e.g., to prevent that the chain (main or side) at side i folds back into itself). To this end, we compute function T ( i, j ) defined as T ( i, j ) = (cid:88) a = { , , , } f a ( i ) f a ( j ) (SI-16)for each pair of beads i and j . T ( i, j ) returns a 1 if and only if the turns t i and t j are along the same axis ( a and ¯ a , respectively). Note that T ( i, j ) is composed of -local terms for the sparse encoding. Firstly, we need to eliminatesequences where the same axis ( a and ¯ a ) is chosen twice in a row (e.g. ) since this will give rise to a chain folding backinto itself. To this end we apply the following penalty term H gc = N − (cid:88) i =3 λ back T ( i, i + 1) (SI-17)with large positive λ back . Note that we can easily control the number of 2-local terms appearing in the sum (linear numberof two local terms). In the general case, for a degree of branching s > , similar terms need to be added to prevent theoverlap of two consecutive bonds within the side chains. t i t i − t i t i − ¯0 ¯1 ¯2 ¯30 ¯3 ¯1 ¯21 ¯2 ¯3 ¯02 ¯3 ¯0 ¯13 ¯1 ¯2 ¯0 TABLE S1. Value of the turn parameter t i (1) imposed by chirality constraint as a function of the values t i − and t i for the case i even (i.e., belonging the the sub-lattice A ). The case i odd is identified by transposing the matrix. Chirality constraints
Natural polymers have a well defined chirality that has to be imposed in our model. In proteins,the position of the side chains at the insertion point with the main chain determines the chirality of each residue. Theposition of the first side chain bead i (1) on the main bead i is imposed by the choice of the (main chain) turns t i − and t i .To enforce the correct chirality, we add a constraint H ch to the Hamiltonian. The required (expected) chirality at i (1) is encoded in the function f ex a ( i (1) ) which is a function of f a ( i − and f a ( i ) defined in Eq.SI-6.Analyzing all possible cases, we can construct the truth-table given in Table S1, which uniquely define f ex a ( i (1) ) for a ∈ { , , , } . To this end, we first need to define an indicator for the parity, g i , which is 0 if i is even and 1 otherwise.We then write g i = − ( − i . Using the information in Table S1 and the definition of g i , we obtain f ex ( i (1) ) = (1 − g i ) (cid:16) f ( i − f ( i ) + f ( i − f ( i ) + f ( i − f ( i ) (cid:17) (SI-18) + g i (cid:16) f ( i ) f ( i −
1) + f ( i ) f ( i −
1) + f ( i ) f ( i − (cid:17) (SI-19) f ex ( i (1) ) = (1 − g i ) (cid:16) f ( i − f ( i ) + f ( i − f ( i ) + f ( i − f ( i ) (cid:17) (SI-20) + g i (cid:16) f ( i ) f ( i −
1) + f ( i ) f ( i −
1) + f ( i ) f ( i − (cid:17) (SI-21) f ex ( i (1) ) = (1 − g i ) (cid:16) f ( i − f ( i ) + f ( i − f ( i ) + f ( i − f ( i ) (cid:17) (SI-22) + g i (cid:16) f ( i ) f ( i −
1) + f ( i ) f ( i −
1) + f ( i ) f ( i − (cid:17) (SI-23) f ex ( i (1) ) = (1 − g i ) (cid:16) f ( i − f ( i ) + f ( i − f ( i ) + f ( i − f ( i ) (cid:17) (SI-24) + g i (cid:16) f ( i ) f ( i −
1) + f ( i ) f ( i −
1) + f ( i ) f ( i − (cid:17) (SI-25)The penalty term to impose the right chirality, H ch , will then the penalize the structures for which f a ( i (1) ) differs from f ex a ( i (1) ) , H ch = λ chirality (cid:88) a N − (cid:88) i =2 (cid:18) − f a ( i (1) ) (cid:19) f ex a ( i (1) ) , (SI-26)for large and positive values of λ chirality . Note that this term actually only contains linear number of 3-local terms for thedenser encoding (5-local for the sparser encoding). C. Construction of H in Definitions.
We can decompose the interaction Hamiltonian into different Hamiltonian terms H ( l ) corresponding to l t h nearest neighbour interactions, H in = H (1) in + H (2) in + H (3) in + . . . (SI-27)To keep the number of terms in H in as low as possible we will need to exploit properties of the tetrahedral lattice. Inthis section, we only consider interactions between beads i and j with j > i + 2 . Indeed, interactions between beads thatare nearest or second nearest neighbour in the primary sequence are present in every conformations so considering themwill not help discriminating the folds. Note also that beads that are separated by less than 5 bonds cannot be nearestneighbour (1-NN) on the lattice. These beads can indifferently represent main chain beads as well as side chain beadssince only the absolute distance measure in bond units along the polymer chain matters (e.g. j (1) − i ≡ j − i + 1 or0 j (1) − i (1) ≡ j − i + 2 ). First, we observe that the parity of j − i corresponds to the parity of d ( i, j ) . Then, since wealternatively choose among positive and negative directions to encode a conformation, we have ( j − i ) = (cid:40) (cid:80) a ∆ n a ( i, j ) = 01 mod 2 if (cid:80) a ∆ n a ( i, j ) = ( − i . (SI-28)In addition, we introduce a notation for the set of first nearest neighbor beads to a given bead i , which will be denoted N ( i ) . For example, N ( i ) = { i − , i + 1 , i (1) } if i is not at the end of the polymer chain or N ( i (1) ) = { i } when limitedto one side chain bead per monomer. Finally, we define a new set of two-index quantities q ( l ) ij , q ( l ) i (1) j , q ( l ) i,j (1) , q ( l ) i (1) j (1) thatencode the presence of an interaction of order l between beads i and j (i.e, i and j are l-NN). First nearest neighbour interactions: H (1) . The quantity q (1) ij describes a contact between a pair of beads i and j that belong to different lattices and that can be any far apart along the polymer chain. If this contact occurs, the energywill be modified by the contact energy (cid:15) (1) ij . The corresponding term in the Hamiltonian is then q (1) ij (cid:15) (1) ij . However, sincewe keep the configuration and the interaction qubit registers independent from each other, before assigning the energycontribution (cid:15) (1) ij we need to make sure that the two beads are indeed at a distance of . The correct energy term istherefore q ( l ) ij = 1 and reads q (1) ij λ ( d ( i, j ) − , which implies that we apply either the stabilization energy (cid:15) (1) ij when both q (1) ij = 1 and d ( i, j ) = 1 , or a penalty contribution weighted by the positive factor λ . Here d ( i, j ) cannot be equal to 0because the two beads must be on different lattices. Note that if r ∈ N ( j ) and i and j are first nearest neighbours, r can either overlap with i or be at a distance of from i . In order to avoid the overlap of different beads, we penalize theoccurrence of the first option by applying a penalty function. This constraint can be written as q (1) ij λ (2 − d ( i, j )) with λ being a large positive value; it will prevent the formation of a local overlaps between polymer beads in the vicinity ofa nearest neighbor contact. Summing over all the possible contacts between backbone beads, we obtain the general formof the first nearest neighbour interaction Hamiltonian H (1) = N − (cid:88) i =1 N (cid:88) j ≥ i +5 j − i =1 mod 2 h (1) i,j (SI-29)with h (1) i,j = q (1) ij (cid:18) (cid:15) (1) ij + λ ( d ( i, j ) −
1) + (cid:88) r ∈N ( j ) λ (2 − d ( i, r )) + (cid:88) m ∈N ( i ) λ (2 − d ( m, j )) (cid:19) (SI-30)The choice of the ratio λ /λ has to be made wisely. In fact, we need that the term within parentheses in Eq. (SI-30) isalways positive when i and j are not in contact. Since |N ( i ) | ≤ and − ( j − i + 1) < d ( i, r ) − d ( i, j ) < ( j − i + 1) when r ∈ N ( j ) , we need λ > j − i + 1) λ + (cid:15) (1) ij . The penalty values λ and λ therefore depend on i and j . Second nearest neighbour interactions: H (2) . In order to take into account second nearest neighbours interactionsbetween beads i and j (which by definition occupy the same sub-lattice), we need to introduce a new set of interactionqubits: q (2) ij . In this case, there are multiple possible configurations that lead to 2-NN interactions, so multiple interactionqubits are needed for each pair ( i, j ) . We deal with these different cases by splitting the Hamiltonian as follows H (2) = H (2) [1] + H (2) [3 ,
3] + H (2) [3 ,
5] + H (2) [5 , . (SI-31) H (2) [1] deals with the case in which a given configuration leads to a direct contact between beads r and i or m and i where ( r, m ) ∈ N ( j ) (see Fig S1), H (2) [1] = N − (cid:88) i =1 N (cid:88) j ≥ i +4 j − i =0 mod 2 h (2) ij [1] (SI-32)where h (2) ij = (cid:88) r ∈N ( j ) q (1) i,r (cid:15) (2) i,j (SI-33)1 (a) (b) (c) (d) FIG. S1. Different contributions to H (2) : (a) and (b) two different configurations contributing to H (2) [1] , (c) a configurationcontributing to H (2) [3 , and (d) a configuration contributing to H (2) [3 , (or H (2) [5 , depending on ordering of the beads r and m ). The remaining configurations can be classified into the following classes depending of the distances between bead i andthe beads r and m , with r, m ∈ N ( j ) (see Fig S1): • [3 ,
3] : d ( i, r ) = 3 and d ( i, m ) = 3 • [3 ,
5] : d ( i, r ) = 3 and d ( i, m ) = 5 • [5 ,
3] : d ( i, r ) = 5 and d ( i, m ) = 3 .We thus introduce three second contact qubits q (2) ij [3 , , q (2) ij [3 , , q (2) ij [5 , and use the same method to build the Hamil-tonian. For instance, H (2) [3 ,
5] = N − (cid:88) i =1 N (cid:88) j ≥ i +4 j − i =0 mod 2 h (2) ij [3 , (SI-34)with h (2) ij [3 ,
5] = (1 − q (1) i,r )(1 − q (1) i,m )q (2) ij [3 , (cid:18) (cid:15) (2) ij + λ ( d ( i, j ) −
2) + λ (3 − d ( i, r )) + λ (5 − d ( i, m )) (cid:19) (SI-35)Note that we also add a check for q (1) i,r (cid:54) = 1 and q (1) i,m (cid:54) = 1 to avoid double counting of the interaction energy (cid:15) (2) i,j . Also inthis case, we choose λ , λ , λ so that the term in the parenthesis is always positive if d ( i, j ) (cid:54) = 2 . Additional qubits arenecessary when dealing with side chains. The nature of these additional terms can be easily elaborated once the gist of themethod is understood. To take into account l -NN interactions with l > we follow the same method of enumerating allthe possible configurations giving rise to the specific interactions. The number of required qubits will evolve exponentiallywith l but still quadratically with N if truncated to a cutoff maximum value. The interaction energy can also depend onthe configuration of interaction.2 Perdomo et al. Perdomo et al. Babbush et al. Babej et al.
This modelModel Hydrophobic-Polar (HP) Coarse-Grained Coarse-Grained Coarse-Grained Coarse-GrainedLattice all all all all TetrahedralTypes ∞ ∞ ∞ ∞ Interactions Nearest Nearest Nearest Nearest l th NearestLocality log N N N l + 2
Qubits N log( N ) N log( N ) N log( N ) N N exp( l ) Scaling N exp( N ) N log ( N ) exp( N ) N Experiment No D-Wave No Rigetti QPU IBM QPUTABLE S2.
Comparison of existing models for 3D protein native structure prediction using quantum algorithms.
The sore point of the existing models is the locality of the generated Hamiltonian and/or the required number of ancillas (qubitresources). In ‘digital’ quantum computers (last two columns), the coupling between more than two qubits can be decomposed intotwo-qubit controlled operations at the cost of an increased circuit depth. On the other hand, the architecture of analog quantumannealers prevents the optimization of non-Ising like Hamiltonian (i.e., with k -body interactions between qubits). In fact, locality-reduction to an isospectral Ising Hamiltonian ( k = 2 ) requires the introduction of additional ancillas qubits and the use of complexmapping schemes (first three columns). Since k -local terms require k − ancillas to be embedded, the resource requirementsfor both quantum approaches (digital and analog) scales exponentially with the polymer size if the locality of the model dependson the polymer length, N , for coarse-grained models. While for the model presented in this work the locality is independent from N this is not the case for most of the proposed Hamiltonians in the literature (Perdomo et al. , Perdomo et al. , and Babej etal. ). The corresponding scalings and the localities are defined in the text. For our model, they are obtained after constructingthe Hamiltonian on Qiskit. The values recorded are obtained without use of any locality reduction schemes.
II. CONTACT MAP
Contact map for the stabilization of two different secondary structure elements: β -sheet (left) versus α -helix (right).Simulations performed with these contact maps reproduce the correct minimal energy structures (see correspondingcartoons). The α -helix is stabilized by a second line of contacts parallel to the main diagonal, while the (anti-parallel) β -sheet forms a series of contacts along the counter-diagonal. D R V Y V H P F H LD 1 1 1 1 1 1 1 0 0R 1 1 1 1 1 1 1 0 0V 1 1 1 1 1 0 1 1 0Y 0 1 1 1 1 1 1 1 1V 0 0 1 1 1 1 1 1 1H 0 0 1 1 1 1 1 1 0P 0 0 1 1 1 1 1 1 1F 1 1 1 1 1 1 1 1 1H 1 1 1 0 0 0 1 1 1L 1 1 1 0 0 0 0 1 1
FIG. S2.
Contact maps for the stabilization of secondary structure elements: β -sheet (left) and α -helix (right). Theupper triangle of the contact map is designed to stabilize a α -helix turn as showed in the structure on the right, which is obtainedfrom a simulation based on the proposed model Hamiltonian. On the other hand, the contacts in the lower triangle give rise to ananti-parallel β -sheet turn (structure to the left, also obtained from a simulation). The signature of these two secondary structuralelements (a series of parallel and anti-parallel contacts with respect to the diagonal of the contact map) can be clearly identified.Color code: the blue color shades depict nearest, second-nearest and third-nearest neighbour interactions. The grey arrows abovethe sequences point to the stabilizing interactions. III. CONVERGENCE OF CVAR VQE FOR NOISY SIMULATION
Simulation of the folding of the peptide Angiotensin using a realistic noise model for hardware.
Generation . . . . . % o f p o pu l a t i o n -9.88-8.18-6.48-6.29-1.7 0 50 Generation . . . . . m a x s p
Evolution of the ( P ( p ) > ) throughout the minimization process for several low-energy conformations. The corresponding conformations can be found in the main text. Results refer to the noisy simulations with 1024 shots (red barsin Figure 2 of the main article). (b) Same for max p P ( p ) with n s = 1024 and n s = 128 . ∗ [email protected] A. Perdomo, C. Truncik, I. Tubert-Brohman, G. Rose, and A. Aspuru-Guzik, Physical Review A (2008). A. Perdomo-Ortiz, N. Dickson, M. Drew-Brook, G. Rose, and A. Aspuru-Guzik, Scientific Reports (2012), 10.1038/srep00571. R. Babbush, A. Perdomo-Ortiz, B. O’Gorman, W. Macready, and A. Aspuru-Guzik, arXiv preprint arXiv:1211.3422 (2012). T. Babej, C. Ing, and M. Fingerhuth, arXiv:1811.00713 [quant-ph] (2018), 00002 arXiv: 1811.00713. M. Fingerhuth, T. Babej, and C. Ing, arXiv:1810.13411 [quant-ph] (2018), 00002 arXiv: 1810.13411. J. Kempe, A. Kitaev, and O. Regev, SIAM J. Comput. , 30 (2006). P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, N. Moll, G. Salis, A. Fuhrer, M. Ganzhorn, D. J. Egger, M. Troyer, A. Mezzacapo,S. Filipp, and I. Tavernelli, Phys. Rev. A , 022322 (2018).8