Using quantum annealing to design lattice proteins

Quantum annealing has shown promise for finding solutions to difficult optimization problems, including protein folding. Recently, we used the D-Wave Advantage quantum annealer to explore the folding problem in a coarse-grained lattice model, the HP model, in which amino acids are classified into two broad groups: hydrophobic (H) and polar (P). Using a set of 22 HP sequences with up to 64 amino acids, we demonstrated the fast and consistent identification of the correct HP model ground states using the D-Wave hybrid quantum-classical solver. An equally relevant biophysical challenge, called the protein design problem, is the inverse of the above, where the task is to predict protein sequences that fold to a given structure. Here, we approach the design problem by a two-step procedure, implemented and executed on a D-Wave machine. In the first step, we perform a pure sequence-space search by varying the type of amino acid at each sequence position, and seek sequences which minimize the HP-model energy of the target structure. After mapping this task onto an Ising spin glass representation, we employ a hybrid quantum-classical solver to deliver energy-optimal sequences for structures with 30-64 amino acids, with a 100% success rate. In the second step, we filter the optimized sequences from the first step according to their ability to fold to the intended structure. In addition, we try solving the sequence optimization problem using only the QPU, which confines us to sizes $\le$20, due to exponentially decreasing success rates. To shed light on the pure QPU results, we investigate the effects of control errors caused by an imperfect implementation of the intended Hamiltonian on the QPU, by numerically analyzing the Schr\"odinger equation. We find that the simulated success rates in the presence of control noise semi-quantitatively reproduce the modest pure QPU results for larger chains.


I. INTRODUCTION
Quantum annealing (QA) [1][2][3][4] is a promising method for finding good solutions to difficult optimization problems.In this method, one aims to find the solution to an optimization problem by encoding it into the ground state of a spin Hamiltonian.The approach of mapping optimization problems to spin systems is not new.It was used already in the 1980s in the context of neural networks [5,6].By exploiting quantum fluctuations and quantum tunneling, QA offers a potentially much faster method for minimizing spin system energies.Technological advances like the D-Wave Advantage quantum annealer, with over 5,000 quantum bits (qubits) and an average connectivity of 15 [7], permit exploration of the QA approach for a wide range of scientifically interesting problems, as illustrated by recent polymer and spin glass studies [8,9].While most of the problems studied concern optimization, attempts are also being made to use quantum annealers for sampling [10][11][12].
We have recently employed this machine for protein folding using the lattice-based HP (hydrophobic/polar) model [13] as a testbed [14].Earlier attempts to use quantum computing for similar folding problems were based upon chain-growth, or turn-based, algorithms [15,16], with which non-local interactions like chain self-avoidance are challenging to implement unless the chain length N is very short (N < ∼ 10).In Ref. [14], we developed a scalable field-like representation, with qubits at all lattice sites, which made it possible to tackle chain lengths up to N = 64 using the D-Wave hybrid quantum-classical solver.
To ensure proper chain configurations, the approach requires penalty terms with Lagrange parameters, but is robust with respect to the choice of these parameters.
Another equally relevant problem from the bio sector is the inverse folding problem [17][18][19], known as protein design, where one seeks to identify a priori unknown sequences which fold into a given structure.Since the structure of a protein is crucial for its function, this problem is highly relevant for drug design.It represents a computational challenge, as it involves the exploration of both sequence and structure spaces.
Here, we tackle the design problem for lattice proteins using QA.We use the HP model as a testbed, for which exact results are available for chain lengths N ≤ 30 [20,21].As is commonly done, we split the design problem into two steps.First, we generate candidate sequences that minimize the energy in the target structure, which we refer to as sequence optimization.Second, we determine whether or not the optimized sequences actually fold into the target structure by either folding the sequence using our aforementioned folding method, or by checking in the databank of all solutions [20,21].Prior work proposed both QA [22] and gate-based [23] schemes for the first of these subproblems, sequence optimization.Here, we present and test a complete QA approach to lattice protein design, comprising both of the above steps.
We find that the hybrid D-Wave annealer can efficiently handle both steps, and thus provides a fast and robust approach to the HP design problem.By contrast, relying entirely on the quantum processing unit (QPU) yields modest performance for larger chains.In order to understand this limitation, we develop and perform time-dependent Schrödinger equation simulations for the sequence optimization problem, using classical high-performance computing clusters.One potential cause of the decrease in success rate with problem size is control errors in the Hamiltonian, which may alter the ground state and thereby lead to incorrect solutions.Indeed, when adding control noise with its strength guided by hardware data, we obtain results that qualitatively reproduce the modest D-Wave pure QPU results.
The approach throughout this work has been evaluated using the D-Wave Advantage as the latter is, currently, the only available quantum annealer with an adequate qubit count.
However, our methods should be valid for any quantum annealer.

II. METHODS
In protein design, one seeks a sequence s = (s 1 , . . ., s N ) that folds into a given target structure, C t .In general, multiple such sequences can exist, and any such sequence constitutes a valid solution.The probability of finding the chain with sequence s in the state C t can be written as where E(C t , s) is the energy of the sequence s in state C t (see Sec. II B), β is inverse temperature and the sum runs over all possible structures C. The design problem therefore translates to finding sequences near the maximum of P β (s).Methods for this task have been developed [24,25].However, maximizing P β (s) involves a generally time-consuming search in both sequence and structure spaces.Therefore, a common approach is to first minimize the energy in the target structure over s, E(C t , s), followed by a filtering step to reject candidate sequences which have a higher probability for a different structure C u .Running folding computations, which determine the most probable structure for a given sequence, is sufficient for this purpose.In this work, we address the design problem by this two-step procedure rather than directly maximizing P β (s).

A. HP lattice proteins
We consider the minimal 2D lattice-based HP model for protein folding [13], in which the protein is represented by a self-avoiding chain of N hydrophobic (H) or polar (P) beads that interact through a pairwise contact potential.A contact between two beads is said to occur if they are nearest neighbors on the lattice but not along the chain.The energy function can be written as E HP = −N HH , where N HH is the number of HH contacts [13].This definition renders the formation of a hydrophobic core energetically favorable.The ground state may be degenerate or unique.For a 2D square lattice, it is known from exhaustive enumerations that about 2% of all HP sequences with N ≤ 30 have a unique ground state [20,21].The availability of exact results for all sequences with N ≤ 30 makes the 2D HP model a useful testbed for novel computational approaches.
Despite their simplicity, coarse-grained HP models are still relevant for the qualitative insights they provide into computationally challenging problems, such as protein folding and design (explored here), liquid-liquid phase separation of intrinsically disordered proteins [26,27], and protein evolution modeling [28,29].

B. HP sequence optimization in QUBO form
Given a target structure, C t , we wish to minimize the energy E HP (C t , s) over sequence, s, using a D-Wave quantum annealer.To this end, the problem must be recast in QUBO form, or, equivalently, into an Ising spin glass format.Furthermore, an auxiliary energy term needs to be included, to control the total number of H beads, N H , in a candidate sequence of length N, since the all-H homopolymer sequence constitutes a trivial solution for unbiased The only information needed about the target structure in order to compute E HP (C t , s) is its connectivity matrix w ij , which indicates whether two arbitrary beads i and j are in contact (w ij = 1) or not (w ij = 0).This holds for any model with pairwise contact interactions, irrespective of the dimensionality and the size of the amino acid alphabet.
When using the HP model, a suitable choice of total energy E(s) to minimize is given by where s i describes whether bead i is of type P (s i = 0) or H (s i = 1).In Eq. ( 2), the first term represents E HP (C t , s), whereas the second term biases the total number of H beads toward a preset value, N H .The balance between the two terms is set by the parameter λ.The 0,1 spins of Eq. ( 2) can be easily transformed into Ising ±1 spins without losing the desired quadratic structure of E. This energy function has a much simpler structure than the corresponding one for the folding problem in Ref. [14], requiring only one Lagrange parameter λ instead of the three in the folding study.
This Lagrange parameter λ must be sufficiently large for the generated sequences to acquire the desired composition, as set by N H . On the other hand, if λ is too large, the energy landscape becomes rugged.Examples of how the efficiency of the hybrid quantumclassical solver varies with λ can be found in Sec.III A. All hybrid production runs were carried out using λ = 2.5.For the pure QPU computations, which are limited to smaller systems, it was possible to use a smaller value, set to λ = 1.1 (Sec.III B).
Minimizing E(s) in Eq. ( 2) can be seen as a graph bisection problem.Unlike a spin glass with nearest-neighbor interactions, it is a fully connected system.Note also that compared to the HP folding problem [14], much fewer spin variables are needed, since only the H/P identity of the beads (whose locations and contacts are fixed for the target structure) needs to be encoded.
For all instances studied in the present paper, it is possible to infer the minimum E HP for a given N H , by inspection of the bead-bead contacts present in the target structure (see Appendices A and B).Hence, it is possible to decide whether or not a proposed solution is correct, even without additional calculations.

C. Hybrid quantum-classical computations
As an alternative to pure QPU computation, the D-Wave Advantage system also offers access to a hybrid quantum-classical solver [30].The hybrid approach uses classical solvers while sending suitable subproblems as queries to the QPU, to speed up the execution and improve the solutions of challenging QUBO problems.Given the limited connectivity (15) within the D-Wave Advantage architecture, the hybrid approach is particularly relevant.
With the hybrid solver, it is possible to tackle problem sizes much larger than with pure QPU computation.
Using the hybrid solver, we performed sequence optimization for three target structures with N = 30, N = 50 and N = 64 (see Fig. 1  In the hybrid approach, the run time needs to be chosen with some care.Indeed, in our previous HP folding study [14], the success rate of the hybrid solver for N > 30 was poor for short run times, while rapidly increasing to values close to 100 % once the run time passed a system size dependent threshold.To determine run times for the sequence optimization problem, we performed a set of preliminary runs for our largest target structure (N = 64), using N H = 42.The hybrid solver consistently returned sequences with the known minimum E HP (Appendix A) for run times ranging from 15 s down to the shortest possible time of 3 s, which is also the default run time for the D-Wave Advantage hybrid solver.Therefore, all the production runs were carried out using this default run time.
To test whether or not the generated sequences actually fold to the desired target structures, we need to perform folding calculations.To this end, we also employ the D-Wave hybrid solver given its demonstrated power for the folding problem [14].Here, for a given optimized sequence s o , the energy E HP (C, s o ) was minimized over chain structure C, using the methods and parameters in Ref. [14] and a 10 × 10 grid.Based on the findings in Ref. [14], the run time was set to 4 s, 120 s and 300 s for N = 30, N = 50 and N = 64, respectively.These run times are larger than the threshold times, above which the success rate was shown to be high [14].

D. Pure QPU computations
The Pegasus topology of the D-Wave Advantage QPU connects each of its qubits to 15 others [7].Problems requiring higher connectivity have to be embedded into the Pegasus graph.This embedding is done by forming "chains" of qubits which act as single qubits.The strength of the coupling between the qubits within a chain is a tunable parameter, called the chain strength.This parameter is typically chosen slightly larger than the minimum chain strength needed to avoid having too many chain breaks (see Sec. III C).
D-Wave offers several so-called samplers for finding embeddings into the QPU topology and performing the QPU computation.We used the DWaveCliqueSampler, designed for dense binary quadratic models [31].It has the property that the chains representing logical qubits share a common length, which facilitates the analysis in Sec.III C. With this method, the number of physical qubits was three times the number of logical qubits in almost all instances studied.For the smallest system (with N = 10), this ratio was two instead of three.All the computations used a chain strength between 2.25 and 4.25 (Appendix B).
The annealing time was set to t f = 2000 µs, its maximum allowed value.The number of output reads per run (annealing cycles), which must be < 10 6 /(t f /µs), was set to 100.

E. Time-dependent Schrödinger equation simulations
As will be seen in Sec.III B, the pure QPU performance deteriorates rapidly with system size.In an attempt to understand this phenomenon, we will perform quantum mechanical simulations with the time-dependent Schrödinger equation.
Similarly, the pure QPU performance was also meager in the folding case [14].With its simple form [Eq. ( 2)], the sequence design problem is better suited for analyzing the shortcomings of the pure QPU performance as compared to the folding case.
We consider an N-qubit system governed by a time-dependent Hamiltonian where H D and H P are the driver and problem Hamiltonians, respectively.On a D-Wave annealer, these two terms take the forms where σ x i and σ z i denote Pauli matrices.For specificity, we will assume a linear annealing schedule given by a(t) = 1 − t/t f and b(t) = t/t f , where t f is the annealing time.We will integrate the Schrödinger equation with this Hamiltonian in time t (with h = 1) using the formalism and algorithm described in Appendix C. The units for time and energy are arbitrary, but related through {unit of time} × {unit of energy} = h.
In addition to the driver Hamiltonian H D in Eq. ( 4), we will also consider the so-called XY -mixer [32], given by which is currently not available on D-Wave's annealers.Transitions generated by this mixer have the property of leaving the constrained sum in Eq. ( 2) unchanged.Hence, if the mixer is used, and the initial state is a uniform superposition of all states satisfying the constraint, it would be possible to minimize E HP at a fixed N H without including the constraint term.

F. Testbed -HP target structures
We seek HP sequences that fold to given target structures with 10-30, 50 and 64 beads.
For N ≤ 30, all HP sequences with unique ground states and the corresponding structures are known from exhaustive enumerations [20,21].This data can be used to decide whether or not a generated sequence actually folds to the desired structure.To evaluate 50-and 64-bead sequences, for which no such data are available, we determine minimum-energy structures by using the D-Wave hybrid solver as described in Ref. [14].

III. RESULTS
Given a target structure (C t ) and a composition (N H ), we wish to find minimum-E HP sequences by minimizing the energy E in Eq. ( 2) on a quantum annealer.In all instances studied, the minimum E HP is known (Appendices A and B), so it is possible to decide whether or not an obtained sequence is a correct solution to the sequence optimization problem.
A sequence that minimizes E HP in C t may, however, have the same or even lower energy in other structures C u = C t .Whether or not this is the case can be checked against existing exact results if N ≤ 30 [20,21], or by performing an energy minimization in the structure space using the hybrid quantum-classical computations as described in Ref. [14].
D-Wave offers solvers based entirely on quantum annealing, as well as a hybrid quantumclassical scheme.In this section, we first try out the hybrid quantum-classical approach with success, even for large chains (Sec.III A).After that, in Sec.III B, we examine the effectiveness of pure QPU calculations, without the classical preprocessing involved in the hybrid approach, using smaller target structures.We observe a rapid decrease in the pure QPU hit rate when increasing the size of the structures.In Sec.III C, we attempt to explain this observation by numerically solving the Schrödinger equation on classical computers.

A. Hybrid quantum-classical computations
Using the hybrid solver, we conducted sequence design for the three target structures shown in Fig. 1 with N = 30, 50 and 64, which will be referred to as T 30 , T 50 and T 64 , respectively.For each target structure, we used a few different compositions, N H (Appendix A).
With no exceptions, the hybrid solver generated sequences with the known minimum E HP (Appendix A) in the target structure, with a 100 % success rate.
Of note, there are sequences that minimize E HP in the target structure without folding to this structure.For such a sequence, the target structure may be one of multiple structures in a degenerate ground state.Alternatively, there exists at least one other structure in which E HP is lower than it is in the target structure, so that a subsequent energy minimization in the conformation space yields a different structure.In such cases, those sequences are not solutions of the design problem for the target structure, and are discarded, even when they are valid solutions for the first step of our two-step approach.This is the price we pay for foregoing an expensive simultaneous search in sequence and structure spaces in favor of a sequence space search as the first step.The sequences emerging from sequence space minimization must be filtered by their ability to fold to the target structure.
Note also that our search for candidate sequences below is not exhaustive; further minimum-energy sequences may exist.Our goal is to find some sequence that folds into the target structure, not all such sequences.
1. Target structure T 30 Our first target structure, T 30 (Fig. 1, left panel), is known from exact results [21] to be the unique ground state of >800 HP sequences.Using the hybrid solver, we minimized E in Eq. ( 2), with λ = 2.5, for this structure for several compositions, 12 ≤ N H ≤ 17.For every N H , 10 hybrid runs all successfully returned sequences with the known minimum E HP .
While most of the thus generated sequences had T 30 as their unique ground state [21], some of them (all with N H = 12, 13 or 16) did not.For each of the latter, a search for possible structures with lower energy was performed, using the hybrid solver (Sec.II C).
No such structure was found, which suggests that the ground states for those sequences are degenerate, and T 30 is one of the structures having the lowest energy.

Target structure T 50
The second target structure, T 50 (Fig. 1, mid panel), comes from a study of a Monte Carlo-based sequence design algorithm [24], which actually optimizes the target population, Eq. ( 1), rather than the energy in this structure.The best sequence found in that study contained 31 H beads [24].
Here, we searched for sequences minimizing the energy E HP of the T 50 structure for N H = 29, 30 and 31, using Eq. ( 2) with λ = 2.5 and the hybrid solver.For N H = 31, the solution to the E HP -minimization problem is unique (Appendix A) and given by the sequence identified in Ref. [24].
As in the T 30 case, for every N H , all 10 hybrid runs successfully gave sequences with the known minimum E HP (Appendix A).For N H = 29 and 30, where the minimum-E HP level is degenerate (Appendix A), the number of distinct sequences obtained from the 10 runs were six and four, respectively.
For each of the 11 distinct optimized sequences, we subsequently minimized E HP over structure by a set of 10 hybrid runs (Sec.II C).Six of the sequences turned out to have T 50 as one of the structures at the lowest lying E HP minimum.For the remaining five sequences, T 50 was the only structure at the global minimum of E HP .Among them, was the sequence with N H = 31 from Ref. [24].The results obtained here support the conclusion that T 50 is the unique ground state of the sequence found in Ref. [24], while at the same time finding several new solutions to the design problem for T 50 .

Target structure T 64
The target structure T 64 (Fig. 1, right panel) has the lowest known energy for an HP sequence with 42 H beads that has been extensively studied [14,33,34].Although not a unique energy minimum, all known structures sharing the same energy (E HP = −42) have a similar shape.The structures differ only within the entirely hydrophobic core, where the chain can be rearranged without altering E HP .
Using Eq. ( 2) with λ = 2.5 and the hybrid solver, we minimized E HP in the target structure T 64 for N H = 36, 40 and 42.As for the previous two target structures, the hybrid solver consistently found sequences with the known minimum E HP (Appendix A) for all N H values.In the N H = 36 and 42 cases, where multiple solutions exist (Appendix A), the returned sequence varied from run to run.
As in the case of T 30 and T 50 , we searched for possible structures with lower E HP using a set of 10 hybrid runs (Sec.II C) per optimized sequence.For every sequence, all runs returned structures with the same E HP as the target structure.The results thus suggest that T 64 is a minimum-E HP structure for all these sequences.However, the minimum is not unique for any of the optimized sequences with N H = 40 or 42.For all these sequences, the core of the T 64 structure is entirely hydrophobic, which makes it possible to rearrange the chain without changing E HP .
Interestingly, the situation appears to be different for one of the optimized N H = 36 sequences, shown in Fig. 2, which has P beads at four core positions.For this sequence, all 10 hybrid folding runs returned the target structure.Changing these four beads to P in the otherwise hydrophobic core appears to lift the degeneracy of the minimum-E HP level.
The above results show that the hybrid quantum-classical method efficiently solves the sequence optimization problem for all the systems studied.All these computations were FIG.
2. An optimized HP sequence that appears to have the target structure T 64 as its unique ground state.The sequence is composed of 36 H (filled) and 28 P (open) beads.For the sequence with 42 H beads studied in Refs.[33,34], this structure is one of several with minimum energy.The degeneracy arises because, for that sequence, the chain can be rearranged in the core of structure without altering the energy.
done with the Lagrange parameter in Eq. ( 2) set to λ = 2.5.To gauge the sensitivity of the success rate to changes in λ, we conducted additional sets of hybrid runs for the three systems (T 30 , T 50 and T 64 with N H = 15, 31 and 36, respectively), using the default run time.We found that λ > 0.25 was necessary for any correct solutions to be found in all three cases.The lower limit on λ is problem-dependent but not larger than 2.0 for any of the systems studied in this paper.As shown in [Fig.3(a)], large values for λ also lead to performance degradation, although, there is a wide window of λ values having a 100 % hit rate, showing that no excessive fine-tuning of λ is required.
The low hit rates at large λ in Fig. 3(a) can be improved at the cost of increasing the run time.show that the N = 30 problem is significantly easier than the other two.
To summarize, here we have designed sequences for three target structures using a two- (b) this widely used approach, a pure sequence space search is performed first, where the types of the amino acids at each position in the target structure are treated as the optimization parameters.An optimized sequence from the first stage is accepted as a solution to the design problem only if a subsequent minimization in the conformation space finds the target structure to be the global minimum for that sequence.We stress that, for both tasks, the hybrid solver gave reliably good results, with robustness with respect to the Lagrange parameter settings.

B. Pure quantum computations
In this section, rather than using the D-Wave hybrid solvers as in the previous section, we explore the ability of pure QPU computations to solve the sequence optimization problem for 10 ≤ N ≤ 20 target structures.We pick, as our target, the most designable structure for a given N denoted by T N , where the designability of a structure is defined as the number of sequences sharing it as their unique ground state.This number is known for all structures with N ≤ 30 from exhaustive enumerations [21].The same databank [21] also provides exact answers to whether or not the generated sequences actually fold to the target structures.Same as panel (a), after restricting the analysis to a filtered dataset without chain breaks (77% of the full dataset) and removing systems with a small energy gap (∆E < 1.0).
For each of the target structures T 10 -T 20 , we performed pure QPU computations for one or a few choices of N H (Appendix B), using the DWaveCliqueSampler.As in the case of the structures in Sec.III A, the minimum E HP , given N H , can be inferred from the contacts present in the target structure, and the solution may be unique or degenerate (Appendix B).
The pure QPU computations recovered all possible solutions to the sequence optimization problem for every (T N , N H ) pair.To quantify the success rate of the pure QPU computations, 10,000 annealing cycles were generated for each combination of T N and N H .The fraction of these yielding a correct solution is referred to as the hit rate as was done in Sec.III A. For this purpose, a correct solution is a solution sequence whose energy matches the previously known minimal E HP for the (T N , N H ) pair.The hit rates obtained in the pure QPU runs can be found in Fig. 4(a).
As noted previously, a sequence that minimizes the energy in the target structure, for a given N H , may not fold to that structure.However, for every (T N , N H ) combination studied here, there is at least one solution to the sequence optimization problem that has the target structure T N as its unique ground state.The precise number of such solutions for different (T N , N H ) pairs can be found in Appendix B.
For the generated sequences that did not have the target structure as its unique ground state, we minimized E HP over structure using the hybrid solver (Sec.II C).In these runs, for almost all the sequences, we found other structures with the same energy as the target structure but none with lower energy, suggesting that the target structure is part of a degenerate ground state.The only exceptions occurred for the target structure T 13 and N H = 6.In this case, five of 18 minimum-E HP sequences attained a lower energy in other structures.Nothing in our procedure precludes the existence of lower-energy structures for a sequence obtained by minimizing the target structure energy.However, had such situations been more common, the recipe used here (sequence space optimization followed by filtering based on folding runs) would not be effective.Note that for T 13 and N H = 6, the existence of lower-energy alternatives is intuitively unsurprising, considering that the target structure has a relatively high energy for the given N H (Appendix B).
Although the pure QPU computations recovered all possible solutions to the sequence optimization problem, the hit rate was strongly problem-dependent [Fig.4(a)].While also depending on N H , the hit rate shows a clear decreasing trend with the problem size N, which limits the range of N which is meaningful to study.
There are several factors that may contribute to the rapid decay of the pure QPU hit rate with problem size as seen in Fig. 4(a).These include (i) finite-t f effects, (ii) chain breaks, (iii) thermal noise and (iv) control errors.Here, we briefly comment on factors (iiii), whereas factor (iv), which relates to the implementation of the couplings J ij and h i in Eq. ( 2), will be discussed in Sec.III C below.
(i) Finite t f .An obvious potential source of error is the use of a finite annealing time t f (≤ 2000 µs).The t f -dependence of the pure QPU hit rate is illustrated in Fig. 5(a) by data obtained for the target structure T 12 and two values of N H (4 and 6).In both cases, the hit rate does increase with t f for small t f .However, it levels off around t f = 400 µs for N H = 4 and already before t f = 100 µs for N H = 6, at values well below unity.This behavior suggests that there must be other, more important error sources than the upper limit on t f .This conclusion is further supported by results from numerically integrating the Schrödinger equation (Sec.II E).Here, we computed the ground-state probabilities, P g , for the different systems at t = t f , for a fixed t f [Fig.5(a)].With our choice of t f = 20 a.u., the P g data roughly match the measured pure QPU hit rates for small N [Fig.4(a)] but not the rapid decay with N seen in the latter case.
(ii) Chain breaks.When embedding the problems into the Pegasus graph, the pure QPU solver creates chains of strongly coupled physical qubits, collectively behaving as a single qubit.In computations, it may happen that such chains, representing logical qubits, break.
To check how our results are affected by such chain breaks, we recomputed the pure QPU hit rates using data only from annealing cycles in which no chain break occurred (77% of the full dataset).A scatter plot comparing the original and recomputed pure QPU hit rates can be found in Appendix D. In many cases, the removal of chain breaks leads to a statistically significant change of the measured hit rate, although the overall agreement between the two datasets is quite good (Pearson correlation coefficient 0.94), (iii) Thermal noise.The effects of thermal noise are likely to be more severe if the energy gap, ∆E, between the ground state and the first excited state of the problem Hamiltonian is small.In our systems, it turns out that ∆E takes on one of three possible values, namely 0.1, 1.0 and 1.1 (Appendix B).In Fig. 4(a) above, there are three systems with markedly lower hit rates than the others, all of which have a small energy gap, ∆E = 0.1.It is conceivable that the low hit rates for these systems, at least in part, is due to thermal noise.However, modeling thermal effects is difficult without extensive details of the underlying physics of the D-Wave processors.
In what follows, we remove from the analysis systems where thermal effects are potentially much stronger than in the others, by focusing on systems with ∆E ≥ 1.0.Furthermore, we will use the filtered dataset without chain breaks, for a cleaner comparison with results from Schrödinger simulations.Redrawing Fig. 4(a) after making these two restrictions, we obtain Fig. 4(b), where the pure QPU hit rate falls off roughly exponentially with N, albeit still with some scatter.
C. Probing the effects of control noise on the pure QPU success rate One potentially limiting factor in the pure QPU computations is analog control errors in the fields h i and couplers J ij [35,36] of the Ising Hamiltonian H P [Eq.( 4)], which are referred to as integrated control errors in D-Wave's documentation [31].The presence of control errors, δh i and δJ ij , leads to a perturbed Hamiltonian whose ground state may not coincide with that of the intended Hamiltonian H P [Eq.( 4)].
Previous work showed that small errors in individual parameters collectively can cause an exponential decay of the success rate with problem size [35,36].
To be able to explore the effects of control errors in our pure QPU computations, we have to make some simplifying assumptions.First, following Refs.[35,36], we assume that all errors δh i and δJ ij are statistically independent and normally distributed, with zero mean and standard deviations σ h and σ J for all δh i and δJ ij , respectively.Second, for computational reasons, we consider only logical qubits, thus essentially ignoring the auxiliary qubits needed when embedding the problems into the QPU topology.However, we take the QPU embedding into account in setting the values of σ h and σ J (see below).
In D-Wave QPU computations, all couplers and fields are rescaled to ĥi = h i /r and Ĵij = J ij /r, where r is the smallest number such that all rescaled parameters fall in given In all systems studied here, the rescaling factor r is set by the chain strength J cs (see Sec. II D) and given by r = J cs /J max .In particular, this implies that the energy gap of the rescaled problem Hamiltonian scales as 1/J cs , which should lead to a decrease in success rate with increasing J cs .
The strength of the control noise on D-Wave's systems has been investigated [31].D-Wave Support suggests using σ h = x max |h i | and σ J = x max |J ij | with x = 0.015, where the maxima are taken over all fields and couplers of the Hamiltonian, including those associated with auxiliary qubits.As indicated above, in our systems, the largest |J ij | is the chain strength J cs .Following the suggestion above, we therefore set σ J = xJ cs .Each logical qubit is represented by a chain of k physical qubits (Sec.II D), where k = 2 for the N = 10 system and k = 3 for all other systems studied.As the physical qubits representing a logical qubit with field h i have fields , where the square root comes from summing over k physical qubits.Summarizing, we then have Using Eq. ( 7) with k and J cs as in the pure QPU computations and x = 0.015, we generated 10,000 perturbed Hamiltonians HP [Eq.( 6)] for each of the systems in Fig. 4(b), and determined the fraction of these, P g , that shared ground state with the intended Hamiltonian H P .For the two systems with (N, N H ) = (10, 4) and (N, N H ) = (20,11), we performed additional sets of pure QPU computations to elucidate how the hit rate depends on J cs .
Figure 6 shows the J cs -dependence of the pure QPU hit rate for these two systems, along with the simulated ground-state probability P g .Clearly, J cs must be sufficiently large to ensure chain stability.On the other hand, we expect the hit rate to drop if J cs gets too large, due to the 1/J cs scaling of the energy gap (see above).The data confirm that J cs must be neither too small nor too large for good performance (Fig. 6), and therefore needs to be chosen with some care.It can also be seen that the hit rate depends only weakly on whether the full dataset or the filtered one without chain breaks is used.As chain breaks do not occur in the simulated systems with only logical qubits, P g does not drop at small J cs .
At large J cs , P g decreases at a rate similar to what is observed for the pure QPU hit rate.more slowly than P g with N. Nevertheless, at a semi-quantitative level, the pure QPU hit rates agree quite well with the P g values obtained using the suggested noise strength x = 0.015.While refraining from attempts to fine-tune x, we note that using x = 0.003 or x = 0.030 leads to, respectively, too high or too low P g values, compared to the QPU hit rates.
In conclusion, the results presented here are consistent with the hypothesis that control errors are an important factor behind the strong N-dependence of the pure QPU hit rates.To fully account for the measured hit rates, and in particular their N H -dependence, additional factors need to be considered, including thermal noise.
There are complementary methods, not explored here, to improve on modest pure QPU hit rates.One is to add a post-processing step, in which the QPU output state is subject to a local optimization with a greedy classical algorithm [37], to drive approximately correct solutions to the desired minimum-energy level.However, in most systems studied here (Fig. 7), all first excited states are local minima left unchanged by this method, which suggests that it is of limited use for our systems.Another approach, called "shimming", is  to refine the calibration of the Hamiltonian by using symmetries of the system [38].These symmetries typically assume zero fields, and coupler values related by a sign flip (±J).These conditions are not met in our systems, which makes useful symmetries hard to find.
Let us finally mention that we also investigated whether the XY -mixer offers an advantage in our problem, by comparing the time evolution of the Schrödinger system when using respectively H D and H XY D in Eqs.(4,5) as drivers.In the cases studied, we found that the improvement is at most marginal.

IV. SUMMARY AND OUTLOOK
Protein design, determining sequences corresponding to a given structure, is a highly relevant biophysical problem.Since both sequence and structure space have to be explored, the problem is very challenging, especially for large chains.We have approached this problem for the HP lattice protein model using quantum annealing, by first minimizing the target structure energy to generate sequences, and then checking if the generated sequences do in fact fold to the target structure.
The approach was evaluated by using the D-Wave Advantage hybrid quantum-classical solver for three structures with chain lengths N = 30, N = 50 and N = 64.Without exceptions, the D-Wave hybrid solver swiftly solves the sequence optimization problem, for which the ground state energy can be deduced (Appendix A).These solutions were then tested for their ability to fold to the target structure, a problem that can also be efficiently tackled using the hybrid solver [14].The two-step procedure was successfully applied to all three target structures.In particular, we identified a previously unknown sequence that appears to have the N = 64 structure as its unique ground state.
In addition, we tested using the D-Wave pure QPU for sequence optimization problems with 10 ≤ N ≤ 20, for which solutions exist in the databank.For all structures, sequences that had the desired structure as their unique ground state were found.However, in line with previous results for the folding problem [14], when using only the QPU for the sequence optimization problem, the performance deteriorates with system size.In order to understand this behavior, which very likely is due to hardware-induced phenomena (noise), we solved Overall, it is clear that the quantum annealer naturally lends itself to protein design.
Using the hybrid quantum-classical annealer, both generating sequences, and subsequently filtering them based on their folding ability, work very well.When using just the QPU, without leveraging D-Wave's hybrid solver, we found that generating sequences with the lowest HP energy for the target structure becomes difficult for larger structures.The results presented lend support to the notion that, at least in part, this difficulty originates from imperfect implementation of the problem Hamiltonian.
Replacing the X-driver with an XY -driver in the Schrödinger simulations, which has been suggested for quantum annealing in general [32], did not improve the results for the sequence optimization problem.This is somewhat surprising since in our case it would remove the only constraint term in Eq. ( 2).
Toward more realistic protein models, one could replace the binary HP encoding spins s i by discrete multi-state spins S i encoding both amino acid type, e.g. in the canonical 20letter alphabet, and the corresponding sidechain conformations, or rotamers [22].Assuming an energy with the quadratic structure E = N i=1 A(S i ) + i<j B(S i , S j ) and a given target backbone conformation, one could then determine the amino acid sequence and rotamers by minimizing E over the S i variables.Given all possible values of all one-and two-body terms A(S i ) and B(S i , S j ), which would have to be predetermined, this minimization could in principle be carried out using QA and a one-hot encoding of the spins S i .However, checking whether or not the optimized sequences fold to the intended and real, rather than lattice-based, backbone conformation would require classical computing.One possibility would be to use the AlphaFold structure prediction method [39,40].Still, due to the Hilbert space dimensionality (2 N ), the numerical implementation requires for a finite-difference ratio, χ(ǫ), which confirm that, for fixed t f , P g displays the expected quadratic dependence on ǫ for small ǫ.

FIG. 1 .
FIG. 1.The 30-, 50-and 64-bead target structures T 30 , T 50 and T 64 used in the sequence design computations with the hybrid solver.

Figure 3 (
b) shows the run time required to attain 50% success rate, τ 1/2 , plotted against λ.For the two larger problems with N = 50 and N = 64, respectively, at λ = 7, the hit rate is tiny when using a run time of 3 s [Fig.3(a)],but can be improved to 50 % by increasing the run time to about 100 s [Fig.3(b)].Note that τ 1/2 grows faster with λ for the two larger systems than it does for the N = 30 system.Not unexpectedly, Figs.3(a,b) both

N H =36 FIG. 3 .
FIG. 3. Dependence of the hit rate on the Lagrange parameter λ [Eq.(2)] when solving the sequence optimization problem by hybrid quantum-classical computations for the target structures T 30 , T 50 and T 64 (Fig. 1) with N H = 15, 31 and 36, respectively.(a) Hit rate as a function of λ when using the default run time, which was 3 s for all systems.Each data point represents an average over 100 runs.(b) The time required to attain 50 % hit rate, τ 1/2 , plotted on a log scale against λ.The horizontal dotted line indicates the default run time (3 s).

FIG. 4 .
FIG. 4. (a) Hit rate when minimizing E(s) in Eq. (2) by pure QPU computation for systems with 10 ≤ N ≤ 20 and different N H (Appendix B).The value of N H is indicated by the plot symbol, and the Lagrange parameter was set to λ = 1.1.For each system, a set of 10,000 annealing cycles was generated using the DWaveCliqueSampler, and the hit rate is the fraction of these that gave a correct solution.The annealing time was set to its maximum value, t f = 2000 µs.The chain strength was chosen individually for each system among the values 2.25, 2.50,. . ., 4.25 for best performance.A circle around the plot symbol indicates that the system has a large energy gap, ∆E ≥ 1.0 (Appendix B).The dashed line is a least-square fit to the encircled data points.(b)

FIG. 5 .
FIG. 5. Effects of using a finite annealing time t f in solving the sequence optimization problem [Eq.(2) with λ = 1.1].(a) Hit rate against t f in pure QPU computations for two of the problems in Fig. 4 (target structure T 12 , N H = 4 and 6).(b) Ground-state probability at t = t f , P g , when numerically integrating the Schrödinger equation (Sec.II E) for systems with 10 ≤ N ≤ 20(Appendix B) using a fixed finite t f (t f = 20 a.u.), plotted against system size, N .The value of N H is indicated by the plot symbol.The Hamiltonian is H(t) = a(t)H D + b(t)H P [Eqs.(3,4)].

Figure 7 (FIG. 6 .
Figure 7(a) shows the simulated hit rates P g , for the systems studied in Fig. 4(b), plotted against N.As in the pure QPU computations [Fig.4(b)], the hit rate falls off roughly exponentially with N. The scatter plot in Fig. 7(b) compares simulated hit rates P g [Fig.7(a)] with pure QPU hit rates [Fig.4(b)].The pure QPU hit rate appears to decrease somewhat

FIG. 7 .
FIG. 7. (a) N -dependence of hit rates from Schrödinger simulations with control noise of strength x = 0.015.The function optimized is E(s) in Eq. (2) with λ = 1.1 and different N H , as indicated by the plot symbol.The systems studied are the same as in Fig. 4(b).For each system, 10,000 perturbed Hamiltonians HP [Eq.(6)] were generated, and the hit rate is the fraction of these that share ground state with the noise-free Hamiltonian, H P .The dashed line is a least-square fit.(b) Scatter plot comparing hit rates from the Schrödinger simulations in panel (a) with the pure QPU hit rates in Fig. 4(b).
the time-dependent Schrödinger equation numerically for different scenarios.Two possible sources of error are inadequate annealing time and control errors in the couplers and fields of the problem Hamiltonian.Whereas inadequate annealing time turned out not to be the problem, our results suggest that control errors have a significant impact on the success rate.Here, we computed the fraction of perturbed Hamiltonians sharing ground state with the original Hamiltonian and compared with the pure QPU results.Employing the same rescaling as on the D-Wave machine and using standard deviations of the noise supplied by D-Wave, we found a semi-quantitative agreement with the observed pure QPU hit rates [Fig.7(b)].A more detailed error model should also incorporate thermal noise.The latter would have to be based on deeper knowledge about the inner workings of the D-Wave architecture, which is, at present, not accessible.These noise investigations were exclusively probing the D-Wave QPU properties for the simple reason that D-Wave is the only available annealer for realistic computations.However, our approach to elucidate the problem by comparing with the time-dependent Schrödinger equation is of wider relevance.

FIG. 11 .
FIG.11.Scatter plot comparing pure QPU hit rates computed using all data [Fig.4(a)] to those obtained after filtering out all annealing cycles in which any chain break occurred.
below).For each target structure C t , we searched for minimum-E HP (C t , s) sequences s, for several fixed compositions, N H .For each combination of C t and N H , we conducted a set of 10 runs, thus generating a set of up to 10 optimized sequences.

TABLE I .
Minimum E HP and the degeneracy of the solution for all sequence optimization instances studied in Sec.III A, using the target structures T 30 , T 50 and T 64 and different compositions, N H .Target structure N H Minimum E HP Degeneracy

TABLE II .
Minimum E HP , the degeneracy of the solution, the energy gap ∆E, and the chain strength J cs used, for all sequence optimization instances studied in Sec.III B, using the target structures T 10 -T 20 and different compositions, N H . Also included is the system used for verification in Appendix C (T 8 ).The number of the solutions that have the target structure as its unique ground state is given within parentheses.∆E is the gap between the two lowest values of E(s)[Eq.(2)], when using λ = 1.1.Target structure N H Minimum E HP Degeneracy ∆E J cs FIG.9.The target structures used for the sequence optimization in Sec.III B and the structure used for verification in Appendix C (T 8 ).The different compositions used can be found in TableII.Appendix C: Numerical integration of the time-dependent Schrödinger equationTo numerically integrate the Schrödinger equation, we first split the time evolution into M steps with length ǫ (t f = Mǫ), yieldingψ(t f ) = U(t M , t M −1 ) ...U(t 1 , t 0 )ψ(0) ,(C1) where ψ(t) is the wave function corresponding to the Hamiltonian in Eq. (3), t m = mǫ (m = 0, . . ., M), and U(t m+1 , t m ) = T exp −i Assuming a linear t-dependence of a(t) and b(t) (Sec.II E), a leading-order Magnus expansion [41] yields ln U(t m+1 , t m ) ≈ −iǫH m + O(ǫ 3 ), where H m = H(t m + ǫ/2).
provides a unitary, second-order accurate integrator (cubic local error).With this approximation, the evolution from time t m to time t m+1 is governed by the constant Hamiltonian H m .