Investigating the potential for a limited quantum speedup on protein lattice problems

Protein folding is a central challenge in computational biology, with important applications in molecular biology, drug discovery and catalyst design. As a hard combinatorial optimisation problem, it has been studied as a potential target problem for quantum annealing. Although several experimental implementations have been discussed in the literature, the computational scaling of these approaches has not been elucidated. In this article, we present a numerical study of quantum annealing applied to a large number of small peptide folding problems, aiming to infer useful insights for near-term applications. We present two conclusions: that even naive quantum annealing, when applied to protein lattice folding, has the potential to outperform classical approaches, and that careful engineering of the Hamiltonians and schedules involved can deliver notable relative improvements for this problem. Overall, our results suggest that quantum algorithms may well offer improvements for problems in the protein folding and structure prediction realm.


Introduction
The structure of a protein captures crucial information about its biological function and therapeutic potential [1]. Knowledge of a proteins' structure unlocks valuable biological information, ranging from the ability to predict protein-protein interactions [2], to structure-based discovery of new drugs [3] and catalysts [4]. Unfortunately, experiments to determine protein structure are challenging and require extensive time and resources [1,5,6]. As of May 2021, the TrEMBL database [7] contained 214 million protein sequences, while only 177,000 protein structures have been deposited in the Protein Data Bank [8]. A reliable computational algorithm for the template-free -that is, when a structurally similar protein, perhaps an evolutionary-related specimen from a different organism, is not available to use as a reference [9] -prediction of a protein's structure and its folding pathway from sequence information alone would enable annotation of millions of proteins and could stimulate major advances in biological research. However, despite steady improvements in the past six decades [10,11,12], and recent advances in the past year [13], a consistent and accurate algorithm for protein folding from sequence has remained elusive.
Over the past decade, there have been attempts to leverage quantum computing for protein structure prediction [14]. The biological structure of a protein is thought to correspond to the minimum of a free energy hypersurface, which for even small peptides is too vast for any classical computer to explore exhaustively [11]. A type of quantum computation that may be appropriate to help is quantum annealing, an approach to exploiting the physics of a controlled quantum system that is considered to be of potential use in optimisation problems (whether classical or quantum in nature) [15]. Typically, the set of possible solutions is mapped to a register of qubits with a binary encoding, and the objective function is represented as a physical Hamiltonian, H problem , whose eigenvectors and eigenvalues are respectively problem solutions and their scores. In particular, the ground state |Φ (or the respective ground eigenspace) corresponds to the global minimum of the problem. The algorithm starts by initialising the register of qubits in the ground state |Ψ(0) of a given Hamiltonian, H trivial , whose ground state is easy to prepare, and gently transforming into the problem Hamiltonian, H problem . If the evolution is slow enough, the adiabatic theorem of quantum mechanics [16] ensures that the final state |Ψ(T ) will be infinitesimally close to |Φ .
Protein chemistry applications of quantum computing have concentrated on a simplified prototype known as the protein lattice model [17], which has been used as a coarse-grained proxy for structure prediction [18,19] (see Figure 1) as well as an invaluable testbed in early theoretical protein folding studies [35]. In this model, the protein is described by a self-avoiding walk on a lattice, whose energy is determined by the contacts between adjacent amino acids, and the minimum energy is identified with the biologically active form of the protein. Unfortunately, the problem of finding the protein configuration that minimises the energy is known to be NP-hard [20,21]. In the context of quantum computing, several encodings (i.e. instructions to map the problem to a Hamiltonian operator and the solutions to a binary string) have been proposed [22,23,24], some of which have also been tested experimentally in D-Wave processors [25,24]. Recent work has attempted extensions of the protein lattice model [26], and even off-lattice models [27,28]. Although multiple algorithmic approaches have been suggested, there is not, to our knowledge, any numerical or analytical study establishing the computational scaling of quantum annealing for protein folding applications.
Protein sequences are constrained to a small range of hundreds, or at most thousands of amino acids. This idiosyncrasy renders considerations of asymptotic Liquorice representation of a randomly generated short protein (peptide) with sequence AVSQVADGILS. In this depiction, every stick represents a bond between two atoms, and the colour of the corresponding half of the stick identifies the nature of the atom: green is carbon, white is hydrogen, blue is nitrogen and red is oxygen. The spheres that surround the sticks, representing van der Waals volume, have been coloured by residue identity. (B) Lattice model of the peptide in (A). The protein is represented as a self-avoiding walk on a lattice, where every node corresponds to a residue. Amino acids that are distant in the sequence but are spatial neighbours induce complex interactions (represented as dotted blue lines) that stabilise a particular fold. Above the dotted lines, we display the Miyazawa-Jernigan stabilisation energy of the contact.
complexity, a typical focus of the quantum algorithms literature, inadequate for a practical evaluation of usefulness. The median length of a human protein is 375 amino acids long [29], while the median clinical target is about 414 amino acids long [30] -and the effective length may be made smaller by considering independently-folding domains, an ubiquitous strategy in computational structure prediction pipelines. Hence, most interesting problems could be encoded with just a few hundred qubits. If there exists a relative, not necessarily asymptotic, speedup with respect to classical heuristics, this will be of interest to both basic research in molecular biology and the pharmaceutical and biotechnological industries, which dedicate significant resources to structural protein studies. Establishing the scaling of quantum annealers in problems of modest size is a crucial baseline for further exploration as near-term quantum annealers capable of addressing larger problems are made available [34].
In this article, we present an extensive numerical study of protein lattice folding in idealised, error-free, closed-system quantum annealing, such as might be achievable in future devices with long coherence times in the presence of error correction [31,32], or with fault-tolerant universal quantum computers employing Hamiltonian simulation [33]. In particular, we have computed the minimum spectral gaps and optimised time-tosolution (TTS) metrics for a large dataset of hard problems. The spectral gaps display a strongly vanishing behaviour, which according to the adiabatic theorem leads one to anticipate a quadratically stronger upper bound in runtime. However, our simulations of unitary evolution reveal a scaling that is several orders of magnitude milder, and that can be optimised further via simple heuristics based on average gap positions.. When compared with classical stochastic optimisation heuristics, we find some evidence of a potential limited quantum speedup for protein lattice folding problems of modest size.

Methodology
Throughout this article, we use a large dataset of peptide (i.e. short protein) sequences, each with a unique global minimum (UGEM), which have been mapped to Ising-like problems with the methods described by Babbush et al. [23]. Peptide sequences with UGEMs have been shown to display the properties of real proteins, even for small sequences (e.g. [35]). The dataset is examined using numerical simulations to assess the gap and expected time-to-solution, studied in terms of the relationship between conditions and results, and compared against an off-the-shelf classical algorithm.

Benchmark dataset generation
We examined a total of 29,503 peptide sequences, an approximately equal number of cases in both dimensions (15,173 in 2D and 14,330 in 3D) and lengths (approximately 4,500 cases per length at a given dimension, with the exception of 2D length 7, where it is challenging to generate UGEM cases, and we considered only 1,700 cases). To produce this dataset, we generated protein sequences by random sampling with replacement of the 20 standard amino acids (ARNDCQEGHILKMFPSTWYV). The states of these instances were enumerated by a brute force algorithm, scoring the energies using the Miyazawa-Jernigan 20-amino acid model [55], and all the sequences with two or more non-equivalent minimum energy conformations were rejected.
These sequences were mapped to an algebraic expression representing the couplings between individual spins in a programmable Ising model. This expression is often known as a Polynomial Unconstrained Binary Optimisation (PUBO) [56]. We employ the turn encoding approach by Babbush et al. [23], which displays the highest reported efficiency in the number of qubits. To perform the mapping in practice, we developed a modified version of SymEngine [57], a computer algebra system (CAS) written in C++, which exploits the idempotency of binary variables leading to up to a five orders of magnitude speedup. We consider protein Hamiltonians with 10 (6 aa peptide in 2D) to 21 qubits (8 aa peptide in 3D), meaning Hilbert spaces of size 1,024 to 2,097,152.

Numerical simulations
2.2.1. Representation of the Hamiltonian Every quantum annealing process considered in this article may be represented in the following general form: where H protein is the Hamiltonian expressing a given protein lattice problem, and H catalyst is one of the following: In these equations, N is the number of qubits; λ is the strength of the catalyst; s is the annealing progress, defined in the interval [0, 1], and generated by some interpolation function s = f (t/T ); I is the identity operator in the N -spin space i.e.
I . In the computational basis, H protein can be represented as a diagonal matrix, and both H start and H stq will be matrices with off-diagonal elements in defined positions. It is easy to prove that the number of null elements in a linear combination of these Hamiltonian matrices grows as O(2 N ), hence significant savings in memory and processing time can be made by exploiting the sparsity of the matrix representation. In contrast, however, H nonstq has a different sparsity pattern and can only be simulated at a much higher computational cost. Our numerical simulations relied on the PETSc library [62,63], where the Hamiltonians were represented in the compressed sparse row (CSR) or Yale format.

Gap evaluation
We estimated the minimum spectral gap between the ground state and the first excited state via numerical diagonalisation of the Hamiltonian matrices, using the Krylov-Schur method [59,60] as implemented in the SLEPc package [61,59]. We computed the eigenvalues in increments of ∆s = 0.1, and interpolated the results using cubic splines. The gap was found as the minimum energy difference between every curve that ran into the ground state at the end of the evolution, and any other curve.
We considered only Hamiltonians without a catalyst term (i.e. λ = 0) for the evaluation of the spectral gap.

Quantum annealing simulation
We estimated the expected time-to-solution (TTS) by simulating the quantum dynamics of the annealer. We assumed an idealised quantum annealer at zero temperature, in the absence of noise, and with perfect control over couplings, which was simulated by numerical integration of the time-dependent Schrödinger equation: where H(t) is the time-dependent Hamiltonian defined in equation 1. We integrated this equation using the Runge-Kutta 5 th order method with adaptive timestepping, as implemented in the PETSc package [62,63]. Runge-Kutta methods have been previously validated for studying quantum annealing [64]. At the end of the evolution, the final state |Ψ(t = T ) is a vector of amplitudes Ψ i whose square modulus |Ψ 2 i | is the probability of measuring a particular binary string in the device.
Several authors have attested the need to use optimal time-to-solution (TTS) metrics to assess the scaling of quantum annealing algorithms, e.g. [43,51]. We employed the Bayesian optimisation package GPyOpt [48] to optimise the annealing time T . We defined an optimisation domain between 0.1 and 1,000 a.u., which was considered acceptable after initial exploration.
We set up a maximum of 50 iterations for annealing programs involving a single parameter (the sample time), and a maximum of 500 iterations for trajectories including a catalyst, where the strength of the catalyst (λ) has to be optimised alongside the sample time; in the case of the non-stoquastic catalyst, given the loss of sparsity of the Hamiltonian matrix, the optimisation was stopped after 30 iterations for some of the length 9 peptides. Default parameters were used otherwise.

Optimal trajectory design
We considered the design of tailored trajectories that increase the efficiency of the algorithm. These trajectories were designed by considering the tendency of gaps to appear more often at certain stages of the annealing trajectory (see Figure S8). We developed a heuristic that accounts for the relative probability of encountering a minimum gap, considering the magnitudes of all the gaps found at that position.
Our heuristic is based in dividing the annealing trajectory in small regions (∆s = 0.05) and estimating the following magnitude for each bin, which we denote "R-score": where k is the index of the bin; the sum extends to all peptides whose minimal gap falls in that bin, indexed by j; P kj is the normalised probability that there is a gap at that position (estimated using kernel density estimation with the Epanechnikov kernel and 0.01 bandwidth); ∆ kj is the magnitude of said gap; and f is a function that weighs the magnitude of the gap into the R-score. The values of the R-score are normalised and used to define a piecewise linear interpolation function that slows down in regions likely to contain a large gap, and anneals at a faster rate otherwise.
We consider two types of optimisation: one employing peptides of all lengths, and a different one designing an optimal program for each length. We also consider several functional forms f (∆ kj ) for the R-score: linear (x), square root ( √ x) and cubic root ( 3 √ x). Performance comparisons for these functions are reported in Figures S9 to S12.

Classical simulated annealing
We compared the performance of the quantum annealing to classical simulated annealing, one of the most common optimisation algorithms in computational biology, which is used in many modern structure prediction packages like Rosetta [49]. We employ an off-the-shelf classical simulated annealing subroutine, the gsl siman.h module of the GNU Scientific Library [50]. This module requires a definition of the optimisation space, and a function that returns the energy of a given configuration. We implemented a simple subroutine to compute the energy of a sequence of lattice moves and, in order to avoid biases deriving from our implementation, we analyse the results in term of Monte Carlo moves instead of times. The classical simulated annealing subroutine employs several parameters, including the number of trials per step, the number of iterations per temperature, the maximum step size of the random walk and the parameters of the Boltzmann distribution (initial and final temperature, Boltzmann constant and damping factor). All of these parameters, except for the Boltzmann constant (set to 1) and the step size (set to 1), were optimised using Bayesian optimisation in the same way parameters were tuned in numerical simulations of quantum annealing.

Statistical analysis of scaling
We assessed the significance of these results using a statistical analysis. We performed a non-linear least squares fit of our data using the lmfit library [65]. We considered four functional models: polynomial (x α ), exponential (e αx ) square exponential (e αx 2 ) and cubic exponential (e αx 3 ). The function was augmented by a constant β, which was set to the average value of the dependent variable for a given length.
We employed three statistical model selection criteria to decide the function that provided the best explanation of the data: the Akaike information criterion (AIC), the Bayesian information criterion (BIC), and the mean squared error (MSE) of the means. The last method was chosen because of the nature of the dataset: the independent variable is discrete, hence the model only provides the prediction for 4 values (3 in 3D). Thus, understanding how the real means differ from the predicted means provides useful information for model selection.

Spectral gap
We start our analysis by studying the minimum spectral gap, ∆, between the ground state and the first excited state. It is often stated (based on theoretical arguments) that the runtime of quantum annealing algorithms is proportional to O(∆ −2 ) [36]. Unfortunately, many problems exhibit an exponentially vanishing gap with increasing problem size [37,38,39], and in particular it is believed that no form of quantum computing is able to efficiently solve NP-complete problems, or at least no such report has withstood scrutiny [40,36]. The distribution of gaps for the protein lattice problem in 2D and 3D, considering a stoquastic process without a catalyst, is shown in Figure  2.
The distribution of gaps at a given length resembles a skewed Gaussian distribution: the majority of gaps are concentrated around a narrow center spanning two orders of magnitude, and there are long, thick tails (containing 5-10% of the data) that spread away several more orders of magnitude to one side. The extent of these tails presents a severe decrease. In the dataset of 2D peptides, the size of the worst-case gaps decreases by five orders of magnitude within a length increase of four amino acids, but only an order of magnitude and a half within the last three lengths considered. Similar results are seen for the 3D peptides. Length 8 three-dimensional peptides show smaller gaps because of their distinct distribution of energetic levels due to cubic symmetry.
This steep decline for the hardest instances is in contrast with the small decrease experienced by an average peptide. The vast majority of the examples populate the area around the median gap, exhibiting a steady but far slower decline. A close examination of the violin plots for the 2D examples also reveals that the position of the peak of the distribution tends to rise as the sequence length increases, and in fact, if we ignore the tails of the distribution, the average gap increases rather than decreases. In interpreting this finding, it is important to keep in mind that our dataset is composed of peptide folding problems that are hard by design, since they have only one ground state solution (plus, in some cases, symmetry-equivalent configurations). These problems are known to be classically very hard [21,20]. In addition to the lack of structure of NP-hard problems, the proportion of lowest-energy solutions in the solution space is minimal, so randomised algorithms will find it very challenging to find the ground state.
We characterised the scaling of the gap using regression analysis by Maximum Likelihood Estimation (see Section 2 for details). Four functional forms were considered: polynomial, x α , exponential, e αx , square exponential, e αx 2 and cubic exponential, e αx 3 . We then employed several standard model selection criteria, detailed in Appendix B, to select the model that better explains the data. The polynomial model x α , with α ≈ −0.75 in 2D and α ≈ −0.4 in 3D, is selected by all criteria, and is significantly better than the second best model.
We also binned the data into different quantiles and repeated the inference, to account for the inhomogeneity of the results (see Tables S1 and S2 in Appendix B). In 2D, the polynomial model is consistently selected across all quantiles, and in 3D, there are some quantiles where the exponential and cubic exponential are selected, which is probably due to the limited range of the dataset and the symmetry effect discussed above. We observed a notable variation of the coefficients across quantiles, as depicted in Figures 2C and 2D. For example, the first quantile, 0-10%, with α = −3.23 in 2D and α = −3.65 in 3D for the polynomial model, contains examples whose gaps vanish at a notably larger rate. On the other hand, the two upper quartiles (75-90% and 90-100%) display positive α values, showing that the gap actually widens with increasing size. Only a portion of the problem instances exhibits fast gap vanishing.
This data does not allow us to conclude that the gap vanishes polynomially. Our results are reminiscent of a previous study on the Exact Cover problem by Young et al. [41,42]. In that study, it was described that, while the scaling of the spectral gap at small problem sizes is consistent with a polynomial [40], at large problem sizes the scaling turns exponential. Similarly, the fraction of problem instances exhibiting small gaps increases at large problem size [42]. We have been unable to study the behaviour of the protein lattice problem at greater sizes, given the large number of qubits and the difficulty of obtaining Hamiltonian expressions beyond 9 amino acids. However, we hypothesise that the protein lattice problem presents exponentially vanishing spectral gaps that will hinder a general polynomial-time solution by quantum annealing for large sizes, while still exhibiting reasonable advantage for problems of modest size, such as may arise in high-throughput peptide discovery experiments.

Simulations of stoquastic dynamics
The minimum spectral gap imposes an upper bound on the running time of quantum annealing, but in order to understand the behaviour of the process we need to access the evolution of the quantum state during the computational process. We investigate this regime by numerical integration of the Schrödinger equation. Since this procedure is costly, the assessment of our entire dataset of circa 30,000 peptides is beyond our resources and, instead, we selected two samples based on spectral gap values. The first sample contains the set of peptides with the smallest gaps (worst-case set), while the second sample is a random selection of peptides (random set). In both cases, each sample contains 100 peptides per chain length of a given dimension, giving a total of 1,400 instances.
A comparison of this sort requires optimising the annealing time to maximise the probability of success. As described by Rønnow et al. [43], a short run can provide a small, but sizeable probability that can be amplified by repetition. In many cases, the repetitions amount to a much shorter runtime than a longer, quasi-adiabatic runtime. We employed Bayesian optimisation [47,48] to find the optimal runtime, as detailed in section 2. The optimised time-to-solution metric, corresponding to the expected runtime to find the correct solution with probability 50%, is shown in Figure 3A-B.
The optimisation of the annealing time of an individual run has a significant impact. We simulated a baseline experiment in which the quantum annealer was run for 1000 a. u., and found an average improvement in the expected total runtime of 15 orders of magnitude in 2D and 10 orders of magnitude in 3D (see Figure S1). We also find a small, but appreciable difference on the dependence on the gap, as depicted in Table 1.
For both the 2D and 3D peptides, the worst-case set containing the smallest gaps does not require significantly longer expected runtimes than the random set. We  Table 1. Correlation coefficients between expected runtime and gap, for the results obtaining optimising the individual runtime and the results obtained with a fixed runtime of 1,000 a.u. ρ is Spearman's rank correlation coefficient between ∆ and T , which describes the monotonic relationship between the two variables (ρ = +1 is perfectly monotonic, ρ = −1 is perfectly inverse monotonic). R is Pearson's correlation coefficient between log 10 1 ∆ 2 and log 10 T performed a two-tailed Welch's t test, and found that the random and worst-case sets could not be distinguished (p-value 0.34, average p-value of subgroup analysis 0.31). In other words, the fact that a problem presents a very small minimum spectral gap does not necessarily indicate it will require a long runtime. This apparently contradictory statement might be explained by the fact that the ∆ −2 scaling in relation to the spectral gap is only valid in the asymptotic limit, and that other terms may govern the dynamics at smaller sizes; and that this rule is above and after all an upper bound, and it is very possible that an instance succeeds even if the annealing is conducted diabatically. In practical terms it is not necessary to guarantee that the final state has a large fidelity with the ground state, but rather to ensure that it has a sizeable amplitude in order to obtain the expected result after enough trials. Another possibility to note is that within the range of problem sizes we inspect there may be only one (or a few) instances of the minimum (or near-minimum) gap occurring during the adiabatic sweep, whereas for large problems a near-minimum gap may occur at multiple points. In addition, the decrease by five orders of magnitude in 2D gaps discussed earlier is also markedly absent from Figures 3A-B We performed a scaling analysis identical to the previous section, finding that, for all cases, either the exponential model is selected over the polynomial, or there is not a significant difference between both of them (see Tables S3 to S6 in Appendix B). In particular, in 2D the polynomial model x α (with α ≈ 0.65) and the exponential model e αx (with α ≈ 0.15 cannot be separated. In 3D, an exponential model e αx with α ≈ 0.45 is selected with high significance. Our findings suggest that the protein lattice problem is not as severely affected by the vanishing of the spectral gap, as might have been expected. Problems with gaps smaller than 10 −2 a.u. (and down to 10 −8 a.u.) do not take significantly longer than problems with a median gap of 0.22 a.u. Moreover, our results hint at an exponential scaling, albeit with a potentially small rate constant. This analysis suggests that this quantum annealing application has a milder scaling than previously expected.

Improvement of the annealing pathway
In the previous section, we found that optimising the sample size can dramatically enhance the performance of quantum annealing, often by several orders of magnitude. This raises the question of whether other simple changes may be able to similarly increase efficiency. In this section we explore several approaches that deliver increasing improvements to the performance of the algorithm.
We first considered whether the linear interpolation function is the most appropriate choice, or if conversely a non-linear interpolation function increases the efficiency of the algorithm. As a first approach, we experimented with three alternative functions: quadratic (x 2 ), cubic (x 3 ) and, in order to include a rapidly changing schedule, sigmoid (1/(1 + e z )). We optimised these functions following the same procedure as in the previous section. The results are summarised in Figures S2 and S3.
Our results suggest that none of the functions considered is able to deliver a significant improvement. The quadratic function displays only minor deviations, and the cubic and sigmoid functions lead to markedly worse performance 99% of the cases. Moreover, when the runtime was reduced, the magnitude of the change was much smaller (roughly 80% of the original time to solution) than when the runtime was slowed down (140-150% in the cubic case and 230-300% in the sigmoid). This supports our initial choice of the linear annealing schedule to establish general trends.
We then considered the introduction of a non-stoquastic catalyst, which has been considered as a potential technique to improve quantum annealing [44,45]. As detailed in the Methods section, this catalyst incorporates terms σ x i σ x j , and is multiplied by a tunable parameter λ that was optimised alongside the sample time. Introducing a catalyst of this form alters the sparsity of the Hamiltonian matrix that our code was optimised for, and therefore we could only consider one of the interpolation functions. This is adequate since our previous examination reveals that the choice of interpolation function has a small impact on efficiency. The results of these simulations are reported in Figures S4 and S5.
The introduction of a non-stoquastic catalyst improves the runtime in about 5% of the worst-case examples, showing a notable average slowdown in the random sample. The magnitude of the deterioration seems to increase with the size of the sequence, while the proportion of sequences that are improved, as well as the magnitude of said improvement, decrease with the size of the sequence. These results suggest that while the non-stoquastic catalyst may well be useful in a fraction of the hardest problems, it is not a general solution. Incidentally, our results agree with those in a recent article by Crosson et al. [46], providing significant analytical and numerical evidence that, in general, stoquastic Hamiltonians are more adept at optimisation than their nonstoquastic counterparts, while allowing however for specific, cleverly engineered nonstoquastic catalysts to provide excellent (sometimes, even exponential) speedups.
The theoretical arguments in [46] led us to consider a stoquastic catalyst that perturbs the annealing trajectory while maintaining stoquasticity. We conceived a catalyst that rotates the transverse field throughout annealing with a sum of terms σ x i . Incidentally, this choice of rotation axis preserves the sparse structure of the matrix and enables rapid simulation. As in the previous case, we introduced a tunable parameter λ that was optimised alongside the sample time.
The stoquastic catalyst also improves the runtime in about 5% of the examples. In this case, however, both the proportion of sequences and the magnitude of the improvement seem to increase with the size of the sequence, and while the improvement is generally more notable in the worst-case dataset than in the random dataset, both present some sort of improvement. The magnitude of the deterioration is significantly smaller than in the non-stoquastic case, and is smaller than an order of magnitude, compared to the 5-6 orders of magnitude introduced before. These results may indicate that a stoquastic catalyst is generally preferable to a non-stoquastic one.
After exploring these alternatives, we considered the design of tailored annealing trajectories: schedules designed to slow down near the spectral gap, and therefore reduce the ratio of excitation to local minima. This is motivated by the tendency of spectral gaps to concentrate around certain particular positions (see Figure S8). In practice, a tailored trajectory is simply a piecewise linear function, defined on equally spaced bins, where the slope of the function at a given bin is inversely proportional to the probability of finding a gap of certain magnitude. We report the heuristic employed to design these schedules in the Methods section. The results of these experiments are summarised in Figures S9 and S10. Tailored trajectories ‡ display better results than any of the previous approaches. Nearly half of the sequences in the random dataset, and about 80% of the worstcase dataset, experience some improvement. The magnitude of the speedup is also significantly larger, with many cases achieving accelerations between 10x and 100x, and we observe an increasing trend where longer peptides improve more than shorter ones. Unfortunately, this technique has a significant caveat: that, when the strategy fails, the deterioration can be far worse (10 4 to 10 6 slow down in some cases).
A primary reason why the tailored approach can result in such slow annealing processes is because our heuristic has to balance the positions where there is a high probability of encountering the gap, and the magnitude of said energy difference. In some cases, particularly in the longer sequences, this balance is failing. A potential solution is to employ the distribution of gaps for every peptide length, instead of considering all the data at once. We therefore designed one annealing schedule per peptide length, and simulated the results as in the previous sections. The results of this experiment are reported in Figures S11 and S12.
This approach is unsurprisingly the most successful of all, improving 60% of the sequences in the random dataset, and an impressive 95% of the sequences in the bad dataset. The magnitude of the improvement is similar, although the average improvement is higher; and we do not see the significant deterioration observed in the previous approach. This results suggest a potential strategy to optimise quantum annealing further.
In practical applications, the position of the gap will not be readily available. However, we have observed that the position of the gaps can be estimated by running short annealing runs with a sudden change at different parts of the schedule. We considered a piecewise linear function, similar to the one used in the tailored trajectories, but changing the slope only at the point where we probe for the gap. When applied to our dataset, this method shows some correlation (full dataset: R = 0.17, p = 10 −6 ; worst-case dataset: R = 0.27, p = 5 · 10 −8 ), suggesting that more involved programs may be able to produce reasonable estimations of the gap position.
The results in this section suggest that simple engineering of the Hamiltonians and trajectories involved in the annealing process can deliver a substantial increase in performance. In this sense, the results presented earlier represent a "sensible baseline" ‡ Here we refer to tailored trajectories employing the cubic root function in the R-score (see Section 2), which are the best performing in our benchmark.
of what is achievable with quantum annealing; and the experiments described in this section propose a pathway to improve the performance further. This also allows us to postulate that further experimentation in this domain may deliver better scalings (in relative terms, but potentially also in complexity) than previously reported.

Comparison with simulated annealing
We have observed that quantum annealing requires exponentially growing runtimes to find the ground state of a protein lattice model, even in the small range of lengths explored in our dataset (see Figure 3). Assuming this remains true of larger systems, as one would expect, it precludes an exponential speedup, since enumerating all possible conformations of a lattice model has the same asymptotic complexity. However, an algorithm which scales significantly better, that is, with a far smaller exponential rate α than the classical case could still be useful for practical applications.
In this section, we compare quantum annealing with classical simulated annealing using the data displayed in Figure 4. Unlike other comparisons of quantum annealing and classical simulated annealing (e.g. [51]), we have not constrained the classical approach to solve a problem in the Ising form. More importantly, we consider a NPhard problem, as opposed to previous work by Albash et al. [51] that considered simpler problems.
The distributions presented in Figure 4A and B show a rapidly growing number of Monte Carlo moves. Visually, the runtime appears to display worse-than-exponential growth. Our model comparison analysis (see Tables S7 to S10 in Appendix B) finds that, in all cases, the model fits to a square exponential e αx 2 with a high level of significance (and this behaviour is reproduced at every quantile).
There are some theoretical arguments (e.g. [52]) which conjecture that simulated annealing converges in exponential time. The square exponential fit found by our statistical analysis could be an artifact of parameter optimisation (note that we optimise four parameters for a single simulated annealing run, as opposed to only one in quantum annealing). This anticipates that quantum annealing provides a better scaling. Our results are made stronger by the fact that in this analysis we have considered only the number of Monte Carlo moves, i.e. the number of evaluations of the energy function. The cost of evaluating this energy is approximately quadratic on the size of the peptide, and the actual performance will be slightly worse than as depicted in Figure 4. Of course, since this cost is polynomial, there will be no changes to the complexity of the algorithm, albeit the practical scaling will be worse.
These findings suggest that quantum annealing has an improved performance over simulated annealing. Over the system sizes we examined, the runtime of quantum annealing scales approximately exponentially, while simulated annealing shows a rapidly growing function that fits better to a square exponential.

Discussion
In this article, we have presented a numerical investigation of quantum annealing applied to protein lattice models. We have considered nearly 30,000 protein sequences, each with a unique global energy minimum, which represent realistic protein problems displaying a folded state, but are also high difficulty instances of a NP-hard problem. We first turned our attention to the minimum spectral gap, a quantity connected theoretically to the runtime of a perfect annealer. We have observed that the gap for these protein sequences can decrease quickly in magnitude, although the scaling appears to be polynomial in the range of sizes considered. The polynomial scaling was confirmed by several statistical selection criteria (detailed in the Methods section and Appendix B), although comparison with prior results reported in the literature leads us to hypothesise that the gap vanishes exponentially. We have also observed that the worst cases decrease by five orders of magnitude between 6 and 9 amino acids. This numerical evidence shows that adiabatic evolution of the computer, where the probability of success nears 100%, will require rapidly growing runtimes and likely be infeasible for worst-case problems.
We then considered optimal annealing runs, where the computer is run for a shorter time to produce a small, but sizeable probability of success that is amplified by repetition. We have established that this runtime grows exponentially with peptide length, although the rate of growth was far smaller than the gap analysis suggested. The scaling was found to be approximately equal to e 0.15L for 2D examples and approximately e 0.75L for 3D examples within the range of problem sizes studied. We also found statistical evidence that peptides with very small gaps are not significantly harder than average cases, and that the exponential rates are almost identical for these two datasets.
We experimented with the parameters of the annealing process in pursuit of strategies to increase efficiency. We considered non-linear interpolation functions, nonstoquastic catalysts, and two adaptive programs to slow down the annealing schedule near regions where a potential gap is expected. Our results demonstrate that the last approach is able to deliver 10-100x speedups. Furthermore, we observe a trend of increasing relative improvements as the size increases: the vast majority of the peptides in the worst-case dataset were improved by this strategy, suggesting that potential decreases in performance can be addressed via careful engineering. We also suggest a method to estimate the gap position when this is not be readily available.
A comparison with classical simulated annealing on our dataset shows that the quantum annealing approach is preferable. Statistical modelling seems to suggest that the scaling of simulated annealing fits best to a square exponential, e αx 2 although theoretical arguments lead us to expect this behaviour to become exponential with a large rate as problem size increases. This implies that for large peptide sizes, a quantum annealer may take significantly less time than a classical machine running a stochastic algorithm.
One of the reasons why quantum annealing may prove useful for protein folding and structure prediction is the limited size of interesting problems. More than half the structures deposited in the Protein DataBank contain fewer than 500 residues, and 80% of the domains in the CATH database [53] are smaller than 200 residues. Similarly, the median length of a human protein is 375 amino acids long [29]. Even if the scaling of quantum annealing is exponential, as long as the exponential rate is low enough to fold small proteins or domains in a timely fashion, this approach will be useful for a multitude of practical problems.
There are other advantages to quantum annealing that could be explored further. When the algorithm fails, the system has been excited to a higher energy state, but although this will only be a local minimum, it may still be useful. For example, if this result is used in a bottom-up approach to explore the conformational space of a protein, it may still be a good starting point for more complex simulations. In contrast, classical simulated annealing is not guaranteed to provide a solution that is close to the global minimum. Future work will investigate this fact, and aim to establish "crossover" points where the scaling worsens.
We believe these results suggest that annealing approaches to quantum optimisation may be a powerful heuristic approach to solve the protein lattice problem, whether in specialised processors [34], near-term digital quantum computers [54] or otherwise. The difficulties of the quantum approach are shared by the classical simulated annealing approach, while the scaling is better. Moreover, even in cases where the annealing approach fails, it can provide solutions that despite not being equivalent to the global minimum, are very close to it, and may be useful in a subsequent refinement procedure. These findings offer encouragement for further research in quantum protein lattice folding and other hybrid quantum-classical algorithms for protein structure prediction.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The code that supports the findings of this study will be made available after acceptance at https://github.com/couteiral/proteinham. Note that some of the results for size 9 peptides had to be restricted to a smaller set of iterations due to increased computational burden, which may explain why size 9 random peptides seem to perform significantly worse than the trend would suggest.