# QFold: Quantum Walks and Deep Learning to Solve Protein Folding

P. A. M. Casares<EMAIL_ADDRESS>Departamento de Física Teórica, Universidad
Complutense de Madrid.  Roberto Campos<EMAIL_ADDRESS>Departamento de Física
Teórica, Universidad Complutense de Madrid.  Quasar Science Resources, SL.  M.
A. Martin-Delgado<EMAIL_ADDRESS>Departamento de Física Teórica, Universidad
Complutense de Madrid.  CCS-Center for Computational Simulation, Universidad
Politécnica de Madrid.

###### Abstract

We develop quantum computational tools to predict how proteins fold in 3D, one
of the most important problems in current biochemical research. We explain how
to combine recent deep learning advances with the well known technique of
quantum walks applied to a Metropolis algorithm. The result, QFold, is a fully
scalable hybrid quantum algorithm that in contrast to previous quantum
approaches does not require a lattice model simplification and instead relies
on the much more realistic assumption of parameterization in terms of torsion
angles of the amino acids. We compare it with its classical analog for
different annealing schedules and find a polynomial quantum advantage, and
validate a proof-of-concept realization of the quantum Metropolis in IBMQ
Casablanca quantum processor.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

Proteins are complex biomolecules, made up of one or several chains of amino
acids, and with a large variety of functions in organisms. Amino acids are 20
compounds made of amine $(-NH_{2})$ and carboxyl $(-COOH)$ groups, with a side
chain that differences them. However, the function of the protein is not only
determined by the amino acid chain, which is relatively simple to figure out
experimentally, but from its spatial folding, which is much more challenging
and expensive to obtain in a laboratory. In fact it is so complicated that the
gap between proteins whose sequence is known and those for which the folding
structure has been additionally analyzed is three orders of magnitude: there
are over 200 million sequences available at the UniProt database
uniprot2019uniprot , but just over 172 thousand whose structure is known, as
given in the Protein Data Bank PDB . Furthermore, experimental techniques
cannot always analyse the tridimensional configuration of the proteins, giving
rise to what is called the dark proteome perdigao2015darkproteome  that
represents a significant fraction of the organisms including humans
bhowmick2016darkproteome ; perdigao2019dark . There are even proteins with
several stable foldings bryan2010metamorphicproteins , and others that have no
stable folding called Intrinsically Disordered dunker2001intrinsically .

Since proteins are such cornerstone biomolecules, and retrieving their folding
so complicated, the problem of protein folding is widely regarded as one of
the most important and hard problems in computational biochemistry, and has
motivated research for decades. Having an efficient and reliable computational
procedure to guess their structure would therefore represent a large boost for
biochemical research.

Until recently, one of the most popular approaches to fold proteins was to
apply a Metropolis algorithm parameterised in terms of the torsion angles, as
is done for example in the popular library package Rosetta Rosetta  and the
distributed computing project Rosetta@Home Rosetta@home ; das2007rosetta@home
. The main problem with this approach, though, is that the problem is
combinatorial in nature, and NP complete even for simple models hart1997robust
; berger1998protein . For this reason, other approaches are also worth
exploring. In the 2018 edition of the Critical Assessment of Techniques for
Protein Structure Prediction (CASP) competition CASP , for example, the winner
was DeepMind’s AlphaFold model AlphaFold , that was able to show that Deep
Learning techniques allow to obtain much better results. DeepMind approach
consisted on training a neural network to produce a mean field potential,
dependent on the distance between amino acids and the torsion angles, that can
be later minimized by gradient descent.

In this article we study how Quantum Computing could help improve the state of
the art in this problem when large error-corrected quantum computers become
available. We propose using the prediction of AlphaFold as a starting point
for a quantum Metropolis-Hastings algorithm. The Metropolis algorithm is a
Markov-chain Monte Carlo algorithm, that is, an algorithm that performs a
random walk $\mathcal{W}$ over a given graph. The Metropolis algorithm is
specially designed to quickly reach the equilibrium state, the state
$\pi^{\beta}$ such that $\mathcal{W}\pi^{\beta}=\pi^{\beta}$. Slowly modifying
the inverse temperature parameter $\beta$ such that the states with smaller
energy become increasingly favoured by the random walk, we should end in the
ground state of the system with high probability.

Figure 1: Example of the smallest dipeptide: the glycylglycine. We can see
that each amino acid has the chain (Nitrogen-$C_{\alpha}$-Carboxy). Different
amino acids would have a different side chain attached to the $C_{\alpha}$
instead of Hydrogen as it is the case for the Glycyne. In each figure we
depict either angle $\phi$ or $\psi$. Angle $\psi$ is defined as the torsion
angle between two planes: the first one defined by the three atoms in the
backbone of the amino acid ($N_{1}$, $C_{\alpha,1}$, $C_{1}$), and the second
by the same atoms except substituting the Nitrogen in that amino acid by the
Nitrogen of the subsequent one: ($C_{\alpha,1}$, $C_{1}$, $N_{2}$). For the
$\phi$ angle the first plane is made out of the three atoms in the amino acid
($N_{2}$, $C_{\alpha,2}$, $C_{2}$) whereas the second plane is defined
substituting the Carboxy atom in the amino acid by the Carboxy from the
preceding amino acid: ($C_{1}$, $N_{2}$, $C_{\alpha,2}$). These graphics were
generated using hanwell2012avogadro  and Inkscape.

Several modifications of the Metropolis-Hastings algorithm to adapt it to a
quantum algorithm have been proposed wocjan2008speedup ; somma2007quantum ;
somma2008quantum ; temme2011quantum ; yung2012quantum ; lemieux2019efficient ,
mostly based on substituting the classical random walk by a Szegedy quantum
walk szegedy2004quantum . On the contrary, our work takes advantage of the
application of a quantum Metropolis algorithm under out-of-equilibrium
conditions similar to what is usually done classically, and has been done on
Ising modelslemieux2019efficient . Specifically, we aim to simulate this
procedure for several small peptides, the smallest proteins with only a few
amino acids; and compare the expected running time with the classical
simulated annealing, and also check whether starting from the initial state
proposed by an algorithm similar to AlphaFold may speed up the simulated
annealing process.

Our work benefits from two different lines of research. The first one makes
use of quantum walks to obtain polynomial quantum advantages, inspired mainly
by Szegedy work szegedy2004quantum , and by theoretical quantum Metropolis
algorithms indicated in the above. In contrast with lemieux2019efficient , our
work focuses only on the unitary heuristic implementation of the Metropolis
algorithm, but studies what happens with a different system (peptides) and
with different annealing schedules instead of only testing a single linear
schedule for the inverse temperature $\beta$. Lastly, we also validate a proof
of concept using IBM Casablanca processor, experimentally realizing the
quantum Metropolis algorithm in actual quantum hardware.

The second line of research related to our work is the use of quantum
techniques to speedup or improve the process of protein folding. The reason
for this is because even simplified models of protein folding are NP hard
combinatorial optimization problems, so polynomial speedups could in principle
be expected from the use of quantum computing instead of their classical
counterparts. The literature on this problem babbush2012construction ;
tavernelli2020resource ; perdomo2012finding ; fingerhuth2018quantum ;
babej2018coarse ; perdomo2008construction ; outeiral2020investigating  and
related ones mulligan2020designing ; banchi2020molecular  focuses on such
simplified lattice models that are still very hard, and mostly on adiabatic
computation. In contrast, our work presents much more realistic fully scalable
model, parametrized in terms of the torsion angles. The torsion angles, also
called dihedral, are angles between the atoms in the backbone structure of the
protein, that determine its folding. An example with the smallest of the
dipeptides, the glycylglycine, can be found in figure 1. These angles are
usually three per amino acid, $\phi$, $\psi$ and $\omega$, but the latter is
almost always fixed at value $\pi$ and for that reason, not commonly taken
into account in the models AlphaFold .

Figure 2: Scheme of the QFold algorithm. Starting from the amino acid
sequence, we use Psi4 to extract the atoms conforming the protein, and a
Minifold module, in substitution of AlphaFold, as initializer. The algorithm
then uses the guessed angles by Minifold as a starting point (or rather, as
the means of the starting von Mises distributions with $\kappa=1$), and the
energy of all possible positions calculated by Psi4, to perform a quantum
Metropolis algorithm that finally outputs the torsion angles. In the scheme of
the algorithm, the backbone builder represents a subroutine that recovers the
links between atoms of the protein, and in particular the backbone chain,
using the atom positions obtained from PubChem using Psi4. The initializer,
instantiated in our case by Minifold, is a second subroutine that gives a
first estimate of the torsion angles, before passing it to the quantum
Metropolis. The energy calculator uses Psi4 to calculate the energy of all
possible rotation angles that we want to explore, and these energies are used
in the quantum Metropolis algorithm, which outputs the expected folding. For a
more detailed flowchart, we refer to the figure 3.

These considerations, and the fact that we use a distilled version of
AlphaFold AlphaFold  as initialization, makes our work different from the
usual approach in quantum protein folding: commonly adiabatic approaches have
been used so far, whereas our algorithm is digital. The downside of this more
precise approach is that the number of amino acids that we are able to
simulate is more restricted, but we are nevertheless able to perform
experiments in actual hardware. Finally, it is worth mentioning that in the
2020 CASP competition, DeepMind’s team tested their AlphaFold v2 algorithm
which significantly improves the results from their previous version.

In summary, the main contributions of our work are threefold: firstly, we
design a quantum algorithm that is scalable and realistic, and provided with a
fault-tolerant quantum computer could become competitive with current state of
the art techniques. Secondly, we analyse the use of different cooling
schedules in out-of-equilibrium quantum walks, and perform ideal quantum
simulations of QFold and compare its performance with the equivalent classical
Metropolis algorithm, pointing towards a quantum speedup. This quantum
advantage is enough to make the quantum Metropolis more convenient than its
classical counterpart in average-length proteins even after taking into
account slowdowns due to error correction protocols. Thirdly, we implement a
proof-of-concept of the quantum Metropolis algorithm in actual quantum
hardware, validating our work.

##  II  QFold algorithm

Figure 3: Flow chart of the QFold algorithm. This figure has to be viewed with
the help of figure 2. QFold has several functionalities integrated altogether,
that could be summarized in an initialization module, a simulation module, and
an experiment module. We denote by diamonds each of the decisions one has to
make. The top part constitutes the initialization module, where Minifold can
be used to get a guess of the correct folding, and Psi4 uses PubChem to
calculate the energies of rotations. The bottom half represents the experiment
or simulation algorithms, that output either Probabilities or
Quantum/Classical TTS and makes use of Qiskit. Different options of the
algorithm are represented by diamonds, and more information on them can be
found in section III.2.

The algorithm we introduce is called QFold and has three main components that
we will introduce in this section (see figure 2 for a scheme of QFold): an
initialization routine to find a good initial guess of the dihedral angles
that characterise the protein folding, a quantum Metropolis to find an even
lower energy state from the initial guess, and a classical metropolis to
compare against the quantum Metropolis to assess possible speedups. The aim of
this section is to introduce the theoretical background we have used for our
results.

###  II.1  Initializer

QFold makes use of quantum walks as a resource to accelerate the exploration
of protein configurations (see figures 2 and 3). However, in nature proteins
do not explore the whole exponentially large space of possible configurations
in order to fold. In a similar fashion, QFold does not aim to explore all
possible configurations, but rather uses a good initialization guess state
based on Deep Learning techniques such as AlphaFold. Since such initial point
is in principle closer to the actual solution in the real space, we expect it
to be most helpful the larger the protein being modelled. In fact, one of the
motivations for our work was the fact that adding a Rosetta relaxation at the
end of the AlphaFold algorithm was able to slightly improve the results of the
AlphaFold algorithm AlphaFold . Notice that a Rosetta relaxation is the way
Rosetta calls its classical Metropolis algorithm. Therefore, we expect that
improved versions of Rosetta, using in our case quantum walks, could be of
help to find an even better solution to the protein folding problem than the
one provided by only using AlphaFold.

The AlphaFold initializer starts from the amino acid sequence ($S$), and
performs the following procedures:

  1. 1.

First perform a Multiple Sequence Alignment (MSA) procedure to extract
features of the protein already observed in other proteins whose folding is
known.

  2. 2.

Then, parametrizing the proteins in terms of their backbone torsion angles
$(\phi,\psi)$ (see figure 1), train a residual convolutional neural network to
predict distances between amino acids, or as they call them, residues.

  3. 3.

Train also a separate model that gives a probability distribution for the
torsion angle conditional on the protein sequence and its previously analysed
MSA features, $P(\phi,\psi|S,MSA(S))$. This is done using a 1-dimensional
pooling layer that takes the predicted distances between amino acids and
outputs different secondary structure such as the $\alpha$-helix or the
$\beta$-sheet 111The $\alpha$-helix and the $\beta$-sheet correspond to two
common structures found in protein folding. Such structures constitute what is
called the secondary structure of the protein, and are characterised because
$(\phi,\psi)=(-\pi/3,-\pi/4)$ in the $\alpha$-helix, and
$(\phi,\psi)=(-3\pi/4,-3\pi/4)$ in the $\beta$-sheet, due to the hydrogen
bonds that happen between backbone amino groups NH and backbone carboxy groups
CO. To make the prediction, the algorithm makes use of bins of size 10º,
effectively discretising its prediction.

  4. 4.

All this information, plus some additional factors extracted from Rosetta, is
used to train an effective potential that aims to give smaller energies to the
configurations that the model believes to be more likely to happen in nature.

Finally, at inference time one starts from a rough guess using the MSA
sequence, and performs gradient descent on the effective potential. One can
also perform several attempts with noisy restarts and return the best option.
Interestingly enough, the neural network is also able to return an estimation
of its uncertainty. Such uncertainty is measured by the parameter $\kappa$ in
the von Mises distribution, and plays the role of the inverse of the variance.
The von Mises distribution is the circular analog of the normal distribution,
and its use is justified because angles are periodic variables
von2014mathematical .

###  II.2  Classical Metropolis

As we have mentioned in the introduction, a relatively popular approach to
perform protein folding has been to use the Metropolis algorithm. The
Metropolis algorithm is an algorithm that performs a random walk over the
configuration space $\Omega$. The configuration space is the abstract space of
possible values the torsion angles that a given protein can take. As such, a
given state $i$ is a list of values for such torsion angles. In particular,
for computational purposes, we will set that angles can take values from a
given set, that is, it will not be a continuous but a discrete distribution.
Over such space we can define probability distributions. Furthermore, since
those angles will dictate the position of the atoms in the protein, the state
$i$ will also imply an energy level $E_{i}$, due to the interaction of the
atoms. In the Rosetta library, the function that calculates an approximation
to such energy is called scoring function.

Starting from a state $i$, the Metropolis algorithm proposes a change
uniformly at random to one of the configurations, $j$, connected to $i$. We
will call $T_{ij}$ to the probability of such proposal. Then this change is
accepted with probability

$A_{ij}=\min\left(1,e^{-\beta(E_{j}-E_{i})}\right),$ (1)

resulting in an overall probability of change $i\rightarrow j$ at a given step
$\mathcal{W}_{ij}=T_{ij}A_{ij}$.

Slowly varying $\beta$ one decreases the probability that steps that increase
the energy of the state are accepted, and as a consequence when $\beta$ is
sufficiently large, the end state is a local minima. If this annealing
procedure is done sufficiently slowly, one can also ensure that the minima
found is the global minima. However, in practice one does not perform this
annealing as slowly as required, resorting instead to heuristic restarts of
the classical walk, and selecting the best result found by the several
different trajectories.

In our implementation we emulate having oracle access to the energies of
different configurations. Such oracle in practice is a subroutine that calls
the Psi4 package turney2012psi4  to calculate the energies of all possible
configurations for the torsion angles that we want to explore. We give more
detail of our particular implementation in section III.2.

###  II.3  Quantum Metropolis

A natural generalisation of the Metropolis algorithm explained in the previous
section is the use of quantum walks instead of random walks. The most popular
quantum walk for this purpose is Szegedy’s szegedy2004quantum , that consists
of two rotations similar to the rotations performed in Grover’s algorithm
Grover . Szegedy’s quantum walk is defined on a bipartite graph. Given the
acceptance probabilities $\mathcal{W}_{ij}=T_{ij}A_{ij}$, $A_{ij}$ defined in
(1), for the transition from state $i$ to state $j$, one defines the unitary

$U\ket{j}\ket{0}:=\ket{j}\sum_{i\in\Omega}\sqrt{\mathcal{W}_{ji}}\ket{i}=\ket{j}\ket{p_{j}}.$
(2)

Taking

$R_{0}:=\mathbf{1}-2\Pi_{0}=\mathbf{1}-2(\mathbf{1}\otimes\ket{0}\bra{0})$ (3)

the reflection over the state $\ket{0}$ in the second subspace, and $S$ the
swap gate that swaps both subspaces, we define the quantum walk step as

$W:=U^{\dagger}SUR_{0}U^{\dagger}SUR_{0}.$ (4)

We refer to Appendix A and Fig. 9 for a detailed account on the use of these
quantum walks in a similar way to the Grover rotations. For completeness, in
Appendix A we review in more detail the theoretical basis of Szegedy quantum
walks that leads us to believe that a quantum advantage is possible in our
problem.

It is well known that if $\delta$ is the eigenvalue gap of the classical walk,
and $\Delta$ the phase gap of the quantum walk, then the complexity of the
classical walk is $O(\delta^{-1})$, the complexity of the quantum algorithm
$O(\Delta^{-1})$, and the relation between the phase and eigenvalue gap is
given by $\Delta=\Omega(\delta^{1/2})$ magniez2011search , offering a
potential quantum advantage. Our algorithm aims to explore what is the
corresponding efficiency gain in practice.

The quantum Metropolis algorithm that we employ lemieux2019efficient  uses a
small modification of the Szegedy quantum walk, substituting the bipartite
graph by a coin. That is, we will have 3 quantum registers: $\ket{\cdot}_{S}$
indicating the current state of the system, $\ket{\cdot}_{M}$ that indexes the
possible moves one may take, and $\ket{\cdot}_{C}$ the coin register. We may
also have ancilla registers $\ket{\cdot}_{A}$. The quantum walk step is then

$\tilde{W}=RV^{\dagger}B^{\dagger}FBV.$ (5)

Here $V$ prepares a superposition over all possible steps one may take in
register $\ket{\cdot}_{M}$, $B$ rotates the coin qubit $\ket{\cdot}_{C}$ to
have amplitude of $\ket{1}_{C}$ corresponding to the acceptance probability
indicated by (1), $F$ changes the $\ket{\cdot}_{S}$ register to the new
configuration (conditioned on the value of $\ket{\cdot}_{M}$ and
$\ket{\cdot}_{C}=\ket{1}_{C}$), and $R$ is a reflection over the state
$\ket{0}_{MCA}$.

Although other clever options are available lemieux2019efficient , here we
implement the simplest heuristic algorithm, which consists of implementing $L$
steps of the quantum walk

$\ket{\psi(L)}:=\tilde{W}_{L}...\tilde{W}_{1}\ket{\pi_{0}},$ (6)

where $t={1,...,L}$ also defines an annealing schedule, for chosen values of
$\beta(t)$ at each step. More detailed explanation of algorithm
lemieux2019efficient  can be found in appendix B.

##  III  Simulations, experiments and results

###  III.1  Figures of merit

When looking for a metric to assess the goodness of given solutions to protein
folding we have to strike a balance between two important aspects: on the one
hand, we want a model that with high probability finds the correct solution.
On the other hand, we would like such procedure to be fast. For example going
through all configuration solutions would be quite accurate albeit extremely
expensive if not directly impossible.

A natural metric to use in this context is then the Total Time to Solution
(TTS) lemieux2019efficient  defined as the average expected time it would take
the procedure to find the solution if we can repeat the procedure in case of
failure:

$TTS(t):=t\frac{\log(1-\delta)}{\log(1-p(t))}.$ (7)

where $t\in\mathbb{N}$ is the number of quantum/random steps performed in an
attempt of the quantum/classical Metropolis algorithm, $p(t)$ the probability
of hitting the right state after those steps in each attempt, and $\delta$ a
target success probability of the algorithm taking into account restarts, that
we set to the usual value of $0.9$. In any case, since it is a constant, the
value of $TTS(t)$ with other value of $\delta$ is straightforward to recover.
Although one should not expect to be able to calculate $p(t)$ in the average
protein because finding the ground state is already very challenging, for
smaller instances it is possible to calculate the $TTS$ for example executing
the algorithm many times and calculating the percentage of them that correctly
identifies the lowest energy state. Using quantum resources, the corresponding
definition is (27) from appendix B.

We can see that this metric represents the compromise between longer walks and
the corresponding expected increase in probability of success. Using this
figure the way we have to compare classical and quantum walks is to compare
the minimum values achieved for the $\min_{t}TTS(t)$. Similar metrics have
also been defined previously in the literature albash2018demonstration .

On the other hand we would also like mention that there is a small
modification of the classical algorithm that could improve its TTS, because we
only output the last configuration of the random walk instead of the state
with minimum energy found so far, a common choice for example in the
Rosetta@Home project. The reason for not having included this modification is
because the length of the classical path, 2 to 50 steps, represents a sizable
portion of the total space that ranges from 64 to 4096 available positions,
whereas that will not be the case for large proteins. We believe that had we
run the classical experiments with that modification, we would have introduced
a relatively large bias in the results, favouring the classical random walks
in the smallest instances of the problem, and therefore likely overestimating
the quantum advantage.

For the experiment run in IBM Quantum systems and whose results can be found
in section III.3.4, the metric we use instead of the TTS is the probability of
measuring the correct answer, and in particular whether we are able to detect
small changes in the probability corresponding to the correct solution.
Measuring the TTS here would not be interesting due to the high level of noise
of the circuit.

###  III.2 Simulation and experimental setup

####  III.2.1  Simulations

For the assessment of our algorithm we have built a simulation pipeline that
allows to perform a variety of options. The main software libraries used are
Psi4 for the calculation of energies of different configurations in peptides
turney2012psi4 , a distilled unofficial version of AlphaFold dubbed Minifold
ericalcaide2019minifold , and Qiskit Qiskit  for the implementation of quantum
walks. Simulations were run on personal laptops and small clusters for
prototyping, and due to its large computational cost, Amazon Web Services aws
for deploying the complete algorithm and obtaining the results. A scheme of
the pipeline of the simulation algorithm can be seen in figure (3).

The initialization procedure takes as input the name of the peptide we want to
simulate and uses Psi4 to download the file of the corresponding molecule from
PubChem kim2019pubchem , an online repository containing abundant information
over many molecules, including atomic positions.

After that, and before executing the quantum and classical walks, the system
checks whether an energy file is available, and if not uses Psi4 to calculate
and store all possible energy values of the different rotations of the torsion
angles. For the calculation of that energy we choose the relatively common
6-31G atomic orbital basis functions jensen2013atomic , and the procedure of
Moller-Plesset to second order helgaker2014molecular  as a more accurate and
not too expensive alternative to the Hartree-Fock procedure. However, it is
also possible to choose any other basis or energy method. Finally, the system
performs the corresponding quantum and random walks and returns the minimum
TTS found.

Figure 4: In section III.3.4 we explain a proof of concept experimental
realization of the quantum Metropolis algorithm, and this figure represents
the implementation of the coin flip subroutine of this hardware-adapted
algorithm, whose operator is represented by $B$ in (5). The circuit has two
key simplifications that reduce the depth. The first one is that we first
rotate the coin register to $\ket{1}$ and then we rotate it back to $\ket{0}$
if the acceptance probability is not 1. This halves the cost, since otherwise
one would have to perform 8 multi-controlled rotations (all possible
combinations of the control values for registers move, $\phi$ and $\psi$), and
in this case we only perform 4 of them. The second simplification again halves
the cost, grouping together rotations with similar values. We empirically see
that the rotation values controlled on $\ket{000}$ and $\ket{010}$ are very
similar, so we group them in rotation $R_{0}$, that implicitly depends on
$\beta$. Similarly, we group the rotations controlled on $\ket{001}$ and
$\ket{101}$ in $R_{1}$, also dependent on $\beta$. This also has the nice
effect of only requiring to control on two out of the three bottom registers,
($\phi$, $\psi$ and move), transforming CCC-$R_{X}$ gates in CC-$R_{X}$.
Separated by the two barriers, from left to right and from top to bottom we
can see the implementation of such CC-$R_{X}$ gates.

In order to evaluate the impact of adding a machine learning module such as
AlphaFold at the beginning of the algorithm, we have implemented the
initialization option to start from random values of the dihedral angles, from
the ones returned by minifold, or the actual original angles that the molecule
has, as returned by the PubChem library.

On the other hand, to evaluate the potential quantum advantage, we also allow
to select the number of bits that specify the discretization of the torsion
angles. For example, 1 bit means that angles can take values in $\\{0,\pi\\}$,
whereas 2 bits indicate discretization in $\pi/2$ radians. In general, the
precision of the angles will be $2^{1-b}\pi$, for $b$ the number of rotation
bits. Notice that the precision of 10º of AlphaFold when reporting their
angles, that we indicated in section II.1, is intermediate between $b=5$ and
$b=6$. The main idea here is that when we increase $b$ or the number of
angles, the size of the search space becomes larger, and evaluating how the
classical and quantum $\min_{t}TTS(t)$ grow we may be able to check whether a
polynomial quantum advantage exists. The $TTS$ can be directly calculated from
inferred classical probabilities, if one is executing the classical metropolis
or the quantum Metropolis in quantum hardware, or from the amplitudes, if one
is running a simulation of the latter. Its specific definition can be seen in
equation (7).

Finally, we implemented and used the experiment mode, that in contrast to the
the simulation mode explained in previous paragraphs, allows us to run the
smallest instances of our problem in actual IBM Q hardware, dipeptides with a
single rotation bit.

Other choices we have to make involve the value of the parameter $\beta$ in
(1), whether it is fixed or follows some annealing schedule, the number of
steps for which the system computes their TTS, or the $\kappa$ parameter from
the von Mises distribution if the initialization is given by Minifold. In
fact, if the original AlphaFold algorithm were to be used, the $\kappa$ values
returned by AlphaFold could actually be used instead of our default value
$\kappa=1$, and furthermore the preparation of the amplitudes corresponding to
this probability distribution could be made efficient using the Grover-Rudolph
state preparation procedure State_prep_grover .

Additionally, while Qiskit allows to recover the probabilities from the
amplitudes, to evaluate the classical walks we have to repeat a certain number
of times the procedure to infer the probabilities. This number of repetitions
is controlled by a variable named number iterations, that we have set to
$500\times(2^{b})^{2n-2}$, where $b$ is the number of bits and $n$ the number
of amino acids, to reflect that larger spaces require more statistics.

Figure 5: Implementation of the full quantum Metropolis circuit from section
III.3.4 implemented in actual quantum hardware, using the coin flip rotation
described in figure 4. The steps of the circuit are separated by barriers for
a better identification (from left to right, and from top to bottom): first
put $\phi$ and $\psi$ in superposition. Then, for each of the two steps
$\tilde{W}$ from (5): put the move register (that controls which angle to
move) in a superposition, operator $V$, and prepare the coin, $B$. Then,
controlled on the coin being in state $\ket{1}$ and the move register the
corresponding angle, change the value of $\psi$ (if move $=\ket{1}$) or $\phi$
(if move $=\ket{0}$), denoted by $F$. Then we uncompute the coin,
$B^{\dagger}$ and the move preparation, $V^{\dagger}$, and perform the phase
flip on $\ket{\text{move}}\ket{\text{coin}}=\ket{00}$ represented by $R$. The
second quantum walk step proceeds equally (with different value of $\beta$ and
the rotations), but now we do not have to uncompute the move and coin
registers before measuring $\phi$ and $\psi$ because it is the last step.

If a non-fixed $\beta$ is attempted, we have implemented and tested several
options for the annealing schedule. The implemented schedules are:

  *   * •

Boltzmann or logarithmic implements the famous logarithmic schedule
kirkpatrick1983optimization

$\beta(t)=\beta(1)\log(te)=\beta(1)\log(t)+\beta(1).$ (8a)

Notice that the multiplication of $t$ times $e$ is necessary in order to make
a fair comparison with the rest of the schedules, so that they all start in
$\beta(1)$. As a consequence, this is not truly the theoretical schedule
required to achieve a quadratic speedup.

  * •

Cauchy or linear implements a schedule given by

$\beta(t)=\beta(1)t.$ (8b)

  * •

geometric defines

$\beta(t)=\beta(1)\alpha^{-t+1},$ (8c)

where $\alpha<1$ is a parameter that we have heuristically set to $0.9$.

  * •

And finally exponential uses

$\beta(t)=\beta(1)\exp(\alpha(t-1)^{1/N}),$ (8d)

where $\alpha$ is again set to $0.9$ and $N$ is the space dimension, which in
this case is equal to the number of torsion angles.

For comparison purposes, the value of $\beta(1)$ chosen has been heuristically
optimized to $50$.

Our current system has two main limitations depending on the mode it is used.
If the aim is to perform a simulation, then the amount of Random Access Memory
of the simulator is the main concern. This is why we have simulated, for a
fixed value of $\beta$:

  * •

Dipeptides with 3 to 5 rotation bits.

  * •

Tripeptides with 2 rotation bits.

  * •

Tetrapeptides with a single rotation bit.

Additionally, dipeptides with 6 bits, tripeptides with 3 bits and
tetrapeptides with 2 bits can be simulated for a few steps, but not enough of
them to calculate our figure of merit with confidence. If the $\beta$ value is
not fixed but follows some annealing schedule, the requirements are larger,
but we can still simulate the same peptides. Further than that, the Random
Access Memory requirements for an ideal (classical) simulation of the quantum
algorithm become quite large. Notice that Qiskit simulator supports 32 qubits
at the moment, but our system is more constrained by the depth of the circuit,
which can run into millions of gates.

####  III.2.2  Experiments

W have also performed experiments in the IBMQ Casablanca processor. In
contrast with the previous experiments, due to the low signal to noise ratio
in the available quantum hardware, it does not make much sense to directly
compare the values of the TTS figure of merit. Instead, the objective here is
to be able to show that we can implement a two-step quantum walk (the minimum
required to produce interference) and still be able to see an increase in
probability associated with the correct state. Since we are heavily
constrained in the depth of the quantum circuit we can implement, we
experiment only with dipeptides, and with 1 bit of precision in the rotation
angles: that is $\phi$ and $\psi$ can be either $0$ or $\pi$.

In this quantum circuit, depicted in figures 4 and 5, we will have 4 qubits,
namely $\ket{\phi}$, $\ket{\psi}$, a coin qubit, and another indicating what
the angle register to update in the next step. Additionally, we always start
from the uniform superposition $\ket{+}_{\phi}\ket{+}_{\psi}$. We then perform
2 quantum walk steps $\tilde{W}$ with values of $\beta$ empirically chosen
$0.1$ and $1$ to have large probabilities of measuring
$\ket{0}_{\phi}\ket{0}_{\psi}$, where we encode the state of minimum energy,
the correct state of the dipeptide.

The figures corresponding to the circuit are 4 and 5, the former depicting the
coin flip procedure and the latter using it as a subroutine in the circuit as
a whole. The coin flip subroutine is the most costly part of the quantum
circuit both in this hardware implementation and in the simulations of
previous sections too, since it includes multiple multi-controlled rotations
of the coin qubit. Perhaps an important remark to make is that this hardware-
adapted circuit contains some simplifications in order to minimize the length
of the circuit as much as possible, since it will be one of the most important
quantities determining the amount of error in the circuit, our limiting
factor.

Peptides | Precision random | Precision minifold | $b$ | quantum min(TTS) random | quantum min(TTS) minifold
---|---|---|---|---|---
Dipeptides | 0.53 | 0.53 | 3 | 136.25 | 270.75
4 | 547.95 | 1137.45
5 | 1426.28 | 1458.02
Tripeptides | 0.46 | 0.71 | 2 | 499.93 | 394.49
Tetrapeptides | 0.51 | 0.79 | 1 | 149.80 | 26.30
Table 1: Table of average precisions defined in equation (9), and
corresponding quantum minimum TTS, defined in equation (7) as the expected
number of steps it would take to find the solution using the quantum
algorithm, with different initializations. $b$ denotes the rotation bits, and
in bold we have indicated which of minifold or random values are best. The aim
of this table is understanding the impact of minifold initialization in the
quantum $\min TTS$, our figure of merit. The two main aspects to notice from
the table are that Minifold precision grows with the size of the peptide, and
that when it is the case that the minifold precision is higher, the
corresponding quantum min TTS values are lower than their random counterparts.
This supports the idea that using a smart initial state helps to find the
native folding of the protein faster.

There is one more important precision to be made: since our implementation of
the quantum circuit in the IBMQ Casablanca processor has 176 basic gates of
depth even after being heavily optimized by Qiskit transpiler, we need a way
to tell whether what we are measuring is only noise or relevant information
survives the noise. Our first attempt to distinguish these two cases was to
use the natural technique of zero-noise extrapolation, where additional gates
are added that do not change the theoretical expected value of the circuit,
but introduce additional noise temme2017error . By measuring how the measured
probabilities change, one can extrapolate backwards to the theoretical ‘zero
noise’ case. Unfortunately, the depth of the circuit is already so large that
it does not work: it does not converge or else returns unrealistic results, at
least when attempted with the software library Mitiq larose2020mitiq .

For this reason we need to find a way out that is only valid because our
circuit is parameterised in terms of the angles and the values of $\beta$.
Additionally we know that if we were to set the value of $\beta$ to $0$, the
theoretical result would be $1/4$ as there are 4 possible states. As a
consequence, our strategy consists of trying to detect changes in the
probability when we use $\beta(\bm{t})=(0,0)$ or $\beta(\bm{t})=(0.1,1)$. The
notation $\beta(\bm{t})$ denotes the value of $\beta$ chosen at each of the
two steps.

For the experiment, we reserved 3 hours of usage of the IBMQ Casablanca
processor, with quantum volume 32. During that time we were able to run 25 so-
called ‘jobs’ with $\beta(\bm{t})=(0,0)$ and 20 ‘jobs’ for 8 arbitrarily
chosen dipeptides. Each run consisted of 8192 repetitions of the circuit
(which can be seen in figures 4 and 5) and an equal number of measurements,
what means that for each dipeptide we run a total of 163840 circuits, and
204800 for $\beta(\bm{t})=(0,0)$ as a baseline.

As we will see from the results, the main limitation of our experiment is the
noise of the system and therefore the depth of the circuit. For this reason we
restrict ourselves to a single rotation bit in dipeptides.

###  III.3  Experimental and simulation results

####  III.3.1 Initialization

In this section we analyse the impact of different initialization methods for
the posterior use of quantum/classical walks. Although we know that AlphaFold
is capable of making a relatively good guess for the correct folding, and
therefore it is reasonable to expect AlphaFold’s guess to be close to the
optimal folding solution in the conformation space, our aim is to give some
additional support to this idea.

As an initializer, we decided to use (and minorly contributed to) Minifold
because even though it does not achieve the state of the art in the prediction
of the angles, it is quite simple and sufficient to illustrate our point.
Minifold uses a residual network implementation given in Tensorflow and Keras
abadi2016tensorflow ; chollet2015keras . Perhaps the most important detail of
using this model is that because we are trying to predict small peptides, and
Minifold uses a window of 34 amino acids for its predictions, we had to use
padding.

The metrics that we analyse in this case are twofold: in the first place, we
would like to see whether Minifold achieves a better precision on the angles
than random guessing. This is a necessary condition for our use of Minifold,
or more generally, any smart initialization module, to make sense. We measure
the precision as $1$ minus the normalized angular distance between the
returned value by the initialization module and the actual value we get from
PubChem (see figures 1 and 2):

$p=1-\frac{d(\alpha,\tilde{\alpha})}{\pi},$ (9)

where $\tilde{\alpha}$ is the estimated angle (either $\phi$ or $\psi$) given
by Minifold or chosen at random, $\alpha$ is the true value calculated from
the output of Psi4 and PubChem, and $d$ denotes the angular distance. Since we
have normalized it, the random initialization gets a theoretical average
precision of 0.5. In table 1 there is a summary comparing the average
precision results of minifold and random initialization broken down by the
protein and bits. The dipeptide results show that due to the small size of the
peptide, minifold has barely a better precision than just random. However,
this situation improves for tripeptides and tetrapeptides, getting a better
precision, and as a consequence, lower TTS values.

Figure 6: Comparison of the Classical and Quantum minimum TTS achieved for the
simulation of the quantum Metropolis algorithm with $\beta=10^{3}$, for 10
dipeptides (with rotation bits $b=3,4,5$), 10 tripeptides ($b=2$) and 4
tetrapeptides ($b=1$), also showing the different initialization options
(random or minifold), and the best fit lines. In dashed grey line we separate
the space where the Quantum TTS is smaller than the Classical TTS. The key
aspect to notice in this graph is that although for smaller instances the
quantum algorithm does not seem to match or beat the times achieved by the
classical Metropolis, due to the exponent being smaller than one (either
$0.89$ or $0.53$ for minifold or random respectively) for average size
proteins we can expect the quantum advantage to be dominant and make the
quantum Metropolis more useful than its classical counterpart. In section
III.3.2 we discuss further details and explain why the random initialization
exponent seems more favourable than the minifold exponent.

Figure 7: Comparison of the Classical and Quantum minimum TTS achieved for the
same peptides as those in figure 6, except that due to computational cost we
do not include dipeptides with 5 rotation bits. This figure corresponds to
section III.3.3, and shows the different initialization options (random or
minifold) and annealing schedules (Boltzmann/logarithmic, Cauchy/linear,
geometric and exponential), and the best fit lines. In dashed grey line we
depict the diagonal. The corresponding fit exponents are given in table 2,
where we can see that in three out of the four cases using an annealing
schedule increases the quantum advantage. On the other hand, using an
exponential schedule does not seem to give but a tiny advantage when used with
a minifold initialization.

If Minifold having a greater precision in the angles than random guessing was
a precondition for our analysis to make sense, the actual metric we are
interested in is whether it has some impact reducing the TTS metric. Otherwise
we could avoid using an initialization module altogether. First of all we were
not expecting an actual reduction of the exponent due to the use of minifold
initialization as much as multiplicative TTS reduction prefactors. However, in
figure 6 the exponent of random initialization model is smaller than the one
corresponding to the minifold initialization in the fit. While this may seem
to indicate that our initialization is in fact harmful to the convergence of
the quantum algorithm, in fact the explanation is quite the opposite: for the
smaller instances of the problem, and very specially in the case of random
initialization, the minimum TTS value is achieved for $t=2$ as can be seen
from the two horizontal structures formed by the blue points in the figure,
meaning that in such cases only using the minifold initialization the quantum
algorithm is able to profit from the incipient quantum advantage. This effect
disappears for larger instances of the problem, but while for random
initialization there is a penalisation of the TTS in the smallest problems
(thus lowering the exponent), the minifold initialization is capable of
correcting for this effect, lowering the TTS of smaller instances of the
problem and, as a bonus, rising the exponent. We are therefore inclined to
believe that the minifold exponent more accurately represents the true
asymptotic exponent of the algorithm for this problem. In conclusion, while
the small size of our experiments does not allow us to see the benefits of
using a smart initialization, they have been important to get a calibrated
estimate of the actual quantum advantage, and we can also see that it helps
reduce the TTS cost both in the classical and quantum algorithm, which nicely
fits the intuition that being closer-than-random to the solution helps to find
the solution faster.

####  III.3.2 Fixed $\beta$

We now discuss whether we are able to observe a quantum advantage in the TTS,
our figure of merit. We will again be discussing the results given in figure
6, for it represents the best fit to the classical vs quantum TTS, and
therefore accurately depicts the expected quantum advantage: we can see that
the slopes separated for the different initialization options are $0.89$ for
the minifold initialization and $0.53$ for the random initialization. As a
consequence, we can see that if these trends are sustained with larger
proteins, there is a polynomial advantage.

As we have seen, figure 6 points towards a quantum advantage. The final
question we would therefore like to answer, is what does this advantage mean
for the modelling of large enough proteins. For that, we only need one
additional ingredient, which is how the expected classical $\min TTS$ scales
with the size of the configuration of the problem. Our data in this respect is
even more restricted because we only have access to configuration spaces of
64, 256 and 1024 positions. Therefore we are making a regression to only three
points, but could give us nevertheless some hints of whether our technique,
the use of quantum Metropolis algorithm, will be helpful to solve the
tridimensional structure of proteins given a large enough and fault tolerant
quantum computer. The regression exponent of a $\log(size)$ vs
$\log(\text{classical }\min TTS)$ fit using both random and minifold
initializations is $r=0.88$, that should not be confused with those in figure
6.

Let us take, for the sake of giving an example, an average 250 amino acids
protein, which has approximately 500 angles to fix. If we use $b=6$ bits to
specify such angles as might be done in a realistic setting, the classical
$\min_{t}TTS$ would be $\approx(2^{b})^{2\times 250\times
r}=(2^{6})^{500\times 0.88}$. The quantum $\min_{t}TTS$, on the other hand,
will be such number to the corresponding exponent of figure 6, that we will
call $e_{m}=0.89$ and $e_{r}=0.53$ for minifold and random. This will
translate to a speedup factor of between $\approx 10^{87}$ and $10^{373}$,
although the latter represents an upper bound and the former is probably
closer to the actual value. This reveals that even if we take into account the
slowdown due to error correction and other factors in quantum computing, the
use of the quantum Metropolis seems very powerful and competitive. And this
conclusion is robust: if we took the quantum advantage exponent to be just
$e=0.95$ and the growth of the TTS with the size of the space an extremely
conservative $r=0.5$, the quantum speedup before factoring in the operating
frequency of the computer will still be a factor of $\approx 10^{22}$. Larger
proteins, which exist in nature, will exhibit even larger speedups in the TTS,
the expected time it would take to find the native structure.

####  III.3.3  Annealing schedules and variable $\beta$

In the previous section we analysed what quantum advantage could be expected
when using a fixed $\beta$ schedule, resulting in a pure quantum walk with
several steps. However, Metropolis algorithms are rarely used in practice with
a fixed $\beta$, since at the beginning of the algorithm one would like to
favour exploration, requiring a low $\beta$ value, and at the end one would
like to focus mostly on exploitation, that is achieved with a high $\beta$
value. It can be seen that this necessity is linked to the well known
exploration-exploitation trade-off in Machine Learning and more generally in
Optimization sutton .

Fit exponents
---
Schedule | Random init. | Minifold init.
Fixed $\beta$ | $0.53$ | $0.89$
Logarithmic | $0.29$ | $0.88$
Linear | $0.34$ | $0.86$
Exponential | $0.37$ | $1.00$
Geometric | $0.29$ | $0.85$
Table 2: Table of scaling exponents for different annealing schedules and
initialization options. The peptides are the same, except that for fixed
$\beta$ we have also included dipeptides with 5 bits of precision, what is
costly for the rest of schedules. For fixed $\beta$, the value heuristically
chosen was $\beta=1000$, while the initial $\beta$ value in each of the
schedules, defined in (8), is $\beta(1)=50$.

The annealing schedule with the strongest theoretical background is usually
called the inhomogeneous algorithm, and it is known because one can prove that
the algorithm will converge to the global minimum value of the energy with
probability 1, albeit generally too slowly (section 3.2 van1987simulated ).
Its implementation is conceptually similar to the Boltzmann or logarithmic
schedule that we use, but with a different prefactor van1987simulated .

Therefore, the question we would like to answer in this section is what
happens when we use our quantum Metropolis algorithm outside of the
equilibrium, that is, when we use a schedule that changes faster than the
theoretical inhomogeneous algorithm. For this task, several optimization
schedules have been proposed, of which we have implemented and tested four
different options whose mathematical formulation can be seen in (8).

Our conclusions from figure 7 and table 2 are that using a variable schedule,
the quantum advantage can be made larger than that of the fixed the
temperature algorithm. However, not all cases give the same advantage, with
the exponential schedule giving practically none, and the geometric schedule
explained in section 5.2 of van1987simulated  being the most promising, with
an exponent of $0.85$ if minifold initialization is used. Linear and
logarithmic schedules lie in between.

Lastly, large differences in the exponents exist depending on the
initialization criteria used, although the same argument on why this is the
case given in section III.3.2 applies here too. In consequence, a more
detailed analysis should be carried out in future work. However, from figures
6 and 7 we can see that the decrease in the expected TTS value when using a
minifold initialization is reflected in the corresponding prefactors, smaller
and more favourable.

####  III.3.4  Experiments in IBMQ Casablanca

Figure 8: Results from hardware measurements corresponding to the experimental
realization of the quantum Metropolis described in section III.3.4 and the
circuit depicted in figures 4 and 5. For each dipeptide we perform a student
t-test to check whether the average success probabilities for
$\beta(\bm{t})=(0,0)$ and $\beta(\bm{t})=(0.1,1)$ are actually different
student1908ttest . The largest p-value measured in all 8 cases is $3.94\cdot
10^{-18}$, indicating that in all cases the difference is significant. For
each of the eight dipeptides we run 163840 times the circuit, and for the
baseline 204800 times.

After having analysed the potential quantum advantage in simulation, we would
also like to implement the quantum Metropolis in IBM Q Casablanca processor.
The circuit we have implemented is made of two quantum walk operators, the
minimum required to cause interference, in the smallest instance of our
problem: dipeptides with a single rotation bit, thus allowing for angles of
$0$ and $\pi$.

The corresponding implementation of the circuit, which was run in IBMQ
Casablanca processor, is depicted in figures 4 and 5. An important detail of
our implementation of this quantum Metropolis algorithm is that in order to be
able to see more than noise we have needed to perform simplifications that
benefit from the particular structure of the problem. This allows us to
minimize the depth of the circuit, that strongly influences the noise, the
limiting figure of our experiment.

In fact, both due to the high level of noise of the circuit, and the minimal
size of the problem we are trying to solve, which only has 4 possible states
to search from, in this experiment it does not make sense to use the Total
Time to Solution as a figure of merit. Rather, we are only interested in
seeing if the quantum circuit is sensitive to the underlying probabilities,
which depend on the chosen $\beta$ value. As we discussed in section III.2.2,
for $\beta=0$ all states are equiprobable whereas for higher $\beta$ values
the probability of the target state is higher, and our experiment aims to see
this difference of probability in practice.

The depicted circuit, with values of $\beta(\bm{t})=(0,0)$ and
$\beta(\bm{t})=(0.1,1)$ was run as explained in the section III.2.2 for
several hundred thousand times for each dipeptide to be able to certify that
our measurements do not correspond to pure noise. The results of such
probability differences are depicted in figure 8, and they are highly
significant as indicated by the low p-value achieved in all cases, smaller or
equal to $3.94\cdot 10^{-18}$, in the corresponding student t-test of the
underlying binomial distribution. Such binomial distribution corresponds to
returning value $1$ when the measurement correctly identifies the minimum-
energy state which is encoded in $\ket{00}$, and $0$ when it does not.

It can also be seen that in seven out of the eight cases tested, the quantum
Metropolis algorithm points in the right direction of increased success
probability. This gives us confidence that we are measuring a small amount of
quantum effects. The outlier, glycylvaline, is surprising because is the
dipeptide that in the simulation shows the greatest theoretical probability of
measuring $\ket{00}$, and at the same time is the dipeptide with the largest
p-value, although still very significant. We can only hypothesise that this is
due to some experimental imperfection.

##  IV  Conclusions and outlook

We have studied how quantum computing might complement modern machine learning
techniques to predict the folding of the proteins. For that, we have
introduced QFold, an algorithm that implements a quantum Metropolis algorithm
using as a starting guess the output of a machine learning algorithm, that in
our case is a simplified implementation of AlphaFold algorithm named Minifold,
but that could in fact be substituted by any initialization module that uses
future improvement of such deep learning techniques.

An important feature of QFold is that it is a realistic description of the
protein folding, meaning that the description of the folded structure relies
on the actual torsion angles that describe the final conformation of the
protein. This is realised by the number of bits $b$ used to fix the precision
of those angles, for which a moderate value of $b=5$ or $b=6$ would be as
accurate as the original AlphaFold. This is in sharp contrast with the rigid
lattice models used to approximate protein folding in the majority of quantum
computing algorithmic proposals for protein folding. Although in our current
simulations presented in this work the range of the precision is limited by
the resources of the classical simulation, nothing prevents QFold from
reaching a realistic accurate precision once a fully fledged quantum computer
is at our disposal, since our formulation is fully scalable within a fault-
tolerant scenario.

The quantum Metropolis algorithm itself relies on the construction of a coined
version of Szegedy quantum walk lemieux2019efficient , and we use the Total
Time to Solution defined in (7), as a figure of merit. Our construction of
this quantum Metropolis algorithm represents the first main contribution of
our work, as it represents a scalable and realistic algorithm for protein
folding that in contrast to previous proposals, does not rely on simplified
lattice models.

The second main contribution is an analysis of the expected quantum advantage
in protein folding, that although moderate in the scaling exponent with the
Minifold initialization, could represent a large difference in the expected
physical time needed to find the minimal energy configuration in proteins of
average size due to the exponential nature of this combinatorial problem. This
quantum advantage analysis is also performed for different realistic annealing
schedules, indicating that the out-of-equilibrium quantum Metropolis
algorithms can show a similar quantum advantage, and can even improve the
advantage shown for the fixed beta case, as can be seen from table 2. The
third contribution is a proof-of-concept small implementation of our algorithm
in actual quantum software.

Our results for the computation of protein folding provide further support to
the development of classical simulations of quantum hardware. A clear message
from our simulations is that it is worthwhile developing quantum software and
classical simulations of ideal quantum computers in order to confirm that
certain quantum algorithms provide a realistic quantum advantage. Some of the
quantum algorithms that must be assessed in realistic conditions are those
based on quantum walks like quantum Metropolis variants with fast annealing
schedule. For them, the complexity scales as $O(\delta^{-a})$, $\delta$ the
eigenvalue gap, and $a>0$ dependent on the specific application and
realization of the quantum Metropolis. This is in contrast with pure quantum
walks, where the classical complexity scales as $O(\delta^{-1})$ and the
quantum complexity as $O(\delta^{-1/2})$, as can be seen in appendix B.
However, it would be naïve to consider this quadratic quantum advantage as an
achievable and useful one in problems similar to ours. Instead, one should aim
to compare quantum algorithms with the classical ones used in practice, where
heuristic annealing schedules are used instead. As a consequence, finding out
the precise values of the corresponding exponent is of great importance since
it is a measure of the quantum advantage. For this reason, not all efforts
should be devoted to finding a quantum advantage solely with noisy ‘baby
quantum computers’, like NISQ devices, but also to continue investigating the
real possibilities of universal quantum computers when they become available.

An unknown that is worth addressing is whether the new software for protein
folding coming out of the CASP competition CASP  will remain in the public
domain or else, will go proprietary. This is specially important since some
the most powerful tools used by that modern software rely on deep learning
that has to be trained. It so happens that the protein databases used for that
training are of public domain. They are the result of many years of
collaborative work among many research institutions that are funded publicly.
Our results point towards the possibility of using an open public software
like Qiskit Qiskit , Psi4 turney2012psi4  and ‘community’ implementations of
the AlphaFold algorithm AlphaFold  of which Minifold ericalcaide2019minifold
is an example, to compensate the power of commercial software for protein
folding.

We would also like to point out that despite the great advances achieved by
the new classical methods in the latest editions of the CASP competition CASP
, there is still huge room for improvement and there are many gaps in the
problem of protein folding that await to be filled. These include
understanding protein-protein and protein-ligand interactions, the real time
folding evolution to the equilibrium configuration, the dark proteome,
proteins with dynamical configurations like the intrinsically disordered
proteins (IDP) and so on and so forth. Crucially, we believe that the current
limitations of the training data sets, which are biased towards easily to be
crystallized proteins, puts constraints on what can be achieved using only
deep learning techniques. Our present research is an attempt to explore
techniques that address these limitations.

Future improvements of our work include primarily a refinement of our
experimentally found quantum advantage, with peptides of larger size, and a
more accurate comparative analysis of the precise quantum advantage that can
be found with each annealing schedule. Such research would be valuable because
it is an important decision to be made when deploying these optimization
algorithms in practice, and no research of this question has been carried out
to the best of our knowledge.

Additionally, we also believe further work should be conducted in clarifying
whether asymptotically one should expect either of the initialization modes to
be polynomially faster finding the ground state. Similar quantum Metropolis
algorithms could be used in a variety of domains, and as such a detailed
analysis, both theoretical and experimental, of the expected quantum advantage
in each case seems desirable.

##  V Acknowledgements

P.A.M.C and R.C contributed equally to this work. We would like to thank kind
advice from Jaime Sevilla on von Mises distributions and statistical t-tests,
Alvaro Martínez del Pozo and Antonio Rey on protein folding, Andrew W. Senior
on minor details of his AlphaFold article, Carmen Recio, Juan Gómez, Juan Cruz
Benito, Kevin Krsulich and Maddy Todd on the usage of Qiskit, and Jessica
Lemieux and the late David Poulin on aspects of the quantum Metropolis
algorithm. We also want to thank IBM Quantum Research for allowing us to use
their quantum processors under the Research program. We also thank Quasar
Science for facilitating the access to the AWS resources. We acknowledge
financial support from the Spanish MINECO grants MINECO/FEDER Projects FIS
2017-91460-EXP, PGC2018-099169-B-I00 FIS-2018 and from CAM/FEDER Project No.
S2018/TCS-4342 (QUITEMAD-CM). The research of M.A.M.-D. has been partially
supported by the U.S. Army Research Office through Grant No. W911NF-14-1-0103.
P. A. M. C. thanks the support of a MECD grant FPU17/03620, and R.C. the
support of a CAM grant IND2019/TIC17146.

##  Appendix A  Szegedy quantum walks

In order to explain what are Quantum Walks, we need to introduce Markov Chain.
Given a configuration space $\Omega$, a Markov Chain is a stochastic model
over $\Omega$, with transition matrix $\mathcal{W}_{ij}$, that specifies the
probability of transition that does not depend on previous states but only on
the present one. Random walks are the process of moving across $\Omega$
according to $\mathcal{W}_{ij}$.

Quantum walks are the quantum equivalent of random walks portugal2013quantum .
The most famous and used quantum walks are those defined by Ambainis
ambainis2007quantum  and Szegedy szegedy2004quantum , although several
posterior generalisations and improvements have been developed such as those
in magniez2011search . Quantum walks often achieve a quadratic advantage in
the hitting time of a target state with respect to the spectral gap, defined
below, and are widely used in several other algorithms paparo2012google ;
paparo2013quantum ; paparo2014quantum , as we shall see.

Szegedy quantum walks are not usually defined using a coin, but rather on a
bipartite walk. This implies duplicating the Hilbert space, and defining the
unitary

$U\ket{j}\ket{0}:=\ket{j}\sum_{i\in\Omega}\sqrt{\mathcal{W}_{ji}}\ket{i}=\ket{j}\ket{p_{j}}$
(10a) and also the closely related
$V\ket{0}\ket{i}:=\sum_{j\in\Omega}\sqrt{\mathcal{W}_{ij}}\ket{j}\ket{i}=\ket{p_{i}}\ket{i}.$
(10b)

Figure 9: Geometrical visualization of a quantum walk operator $W$ of Szegedy
type. $W$ performs a series of rotations that in the subspace
$\mathcal{A}+\mathcal{B}$, defined in (11b) with their corresponding rotation
operators (13b), may be written as a block diagonal matrix, where each block
is a 2-dimensional rotation $\omega_{j}=R(2\varphi_{j})$ given by (18). This
figure represents the direct sum of Grover-like rotations in the subspace
spanned by $\mathcal{A}+\mathcal{B}$, and therefore $W$. This quantum walk
operator represents equation (4) from section II.3.

$\mathcal{W}_{ij}$ must be understood as the probability that state $\ket{i}$
transitions to state $\ket{j}$. We can check that these operators fulfil that
$SU=VS$, for $S$ the Swap operation between the first and second Hilbert
subspace. The cost of applying $U$ is usually called (quantum) update cost.
Define also the subspaces

$\mathcal{A}:=\text{span}\\{\ket{j}\ket{0}:j\in\Omega\\}$ (11a) and
$\mathcal{B}:=U^{\dagger}SU\mathcal{A}=U^{\dagger}VS\mathcal{A}.$ (11b)

Having defined $U$ and $V$, we can define $M:=U^{\dagger}VS$, a matrix with
entries valued
$\braket{i,0}{U^{\dagger}VS}{j,0}=\sqrt{\mathcal{W}_{ji}}\sqrt{\mathcal{W}_{ij}}=\sqrt{\pi_{i}/\pi_{j}}\mathcal{W}_{ij}$,
the last equality thanks to the detailed balance equation 222Notice that in
most texts the definition of $M$ does not explicitly include $S$. It is
assumed implicitly though.. In fact, in matrix terms it is usually written
$M=D_{\pi}^{-1/2}\mathcal{W}D_{\pi}^{1/2}$ where we have written $D_{\pi}$ to
indicate the diagonal matrix with the entries of the equilibrium state-vector
of probabilities, $\pi$. This implies that $\mathcal{W}$ and $M$ have the same
spectrum $\lambda_{0}=1\geq...\geq\lambda_{d-1}\geq 0$, as the matrix
$\mathcal{W}$ is positive definite, $p^{T}\mathcal{W}p\in[0,1]$, and of size
$d$. The corresponding eigenstates are $\ket{\phi_{j}}\ket{0}$, and phases
$\varphi_{j}=\arccos{\lambda_{j}}$. In particular
$\ket{\phi_{0}}=\sum_{i}\sqrt{\pi_{i}}\ket{i}$, is the equilibrium
distribution.

We can also define the projectors around $\mathcal{A}$ and $\mathcal{B}$ as
$\Pi_{\mathcal{A}}$ and $\Pi_{\mathcal{B}}$

$\Pi_{\mathcal{A}}:=(\mathbf{1}\otimes\ket{0}\bra{0}),$ (12a)
$\Pi_{\mathcal{B}}:=U^{\dagger}VS(\mathbf{1}\otimes\ket{0}\bra{0})SV^{\dagger}U$
(12b)

with their corresponding rotations

$R_{\mathcal{A}}=2\Pi_{\mathcal{A}}-\mathbf{1},$ (13a)
$R_{\mathcal{B}}=2\Pi_{\mathcal{B}}-\mathbf{1}.$ (13b)

Using this rotation we further define a quantum walk step as we did in the
main text equation (4),

$W=R_{\mathcal{B}}R_{\mathcal{A}}=U^{\dagger}SUR_{\mathcal{A}}U^{\dagger}SUR_{\mathcal{A}}.$
(14)

Using the previous expressions we can state, using
$\Pi_{\mathcal{A}}\ket{\phi_{j}}\ket{0}=\ket{\phi_{j}}\ket{0}$

$\Pi_{\mathcal{A}}U^{\dagger}VS\ket{\phi_{j}}\ket{0}=\cos\phi_{j}\ket{\phi_{j}}\ket{0},$
(15a) and
$\Pi_{\mathcal{B}}\ket{\phi_{j}}\ket{0}=U^{\dagger}VS\cos\phi_{j}\ket{\phi_{j}}\ket{0}.$
(15b)

(15b) is true because (supplementary material yung2012quantum ),

$\Pi_{\mathcal{A}}U^{\dagger}VS\Pi_{\mathcal{A}}=\Pi_{\mathcal{A}}SV^{\dagger}U\Pi_{\mathcal{A}}$
(16)

due to

$\begin{split}&\braket{\phi_{j},0}{\Pi_{\mathcal{A}}U^{\dagger}VS\Pi_{\mathcal{A}}}{\phi_{j},0}=\lambda_{j}\\\
&=\lambda_{j}^{\dagger}=\braket{\phi_{j},0}{\Pi_{\mathcal{A}}SV^{\dagger}U\Pi_{\mathcal{A}}}{\phi_{j},0}.\end{split}$
(17)

Thus $W$ will preserve the subspace spanned by
$\\{\ket{\phi_{j}}\ket{0},U^{\dagger}VS\ket{\phi_{j}}\ket{0}\\}$, which is
invariant under $\Pi_{\mathcal{A}}$ and $\Pi_{\mathcal{B}}$; mirroring the
situation in the Grover algorithm Grover . Also, as a consequence of the
previous, and of the fact that in $\mathcal{A}+\mathcal{B}$ operator $W$ has
eigenvalues $e^{2i\varphi_{j}}$ magniez2011search ; szegedy2004quantum , in
such subspace the operator $W$ can be written as a block diagonal operator
with matrices

$w_{j}=\begin{pmatrix}\cos(2\varphi_{j})&-\sin(2\varphi_{j})\\\
\sin(2\varphi_{j})&\cos(2\varphi_{j})\end{pmatrix}.$ (18)

Finally, notice that the eigenvalue gap of $\mathcal{W}$ is defined as
$\delta=1-\lambda_{1}$, and in general the hitting time of a classical walk
will grow like $O(\delta^{-1})$ (Proposition 1 of magniez2011search ). On the
other hand, the hitting time of the Quantum Walk will scale like
$O(\Delta^{-1})$ (Theorem 6 of magniez2011search ), where
$\Delta:=2\varphi_{1}$. But $\Delta\geq 2\sqrt{1-|\lambda_{1}|^{2}}\geq
2\sqrt{\delta}$, so $\Delta=\Omega(\delta^{1/2})$. In fact, writing
$\delta=1-\lambda_{1}=1-\cos\varphi_{1}$, and expanding in Taylor series
$\cos\varphi_{1}=1-\frac{\varphi_{1}^{2}}{2}+\frac{\varphi_{1}^{4}}{24}+O(\varphi_{1}^{6})$,
we can see that

$\frac{\varphi_{1}^{2}}{2}\geq
1-\cos\varphi_{1}\geq\frac{\varphi_{1}^{2}}{2}-\frac{\varphi_{1}^{4}}{24}.$
(19)

Using the definitions of $\delta$ and $\Delta$ and the fact that
$\varphi_{1}\in(0,\pi/2)$, it is immediate

$\begin{split}\frac{\Delta^{2}}{8}\geq\delta\geq\frac{\Delta^{2}}{8}-\frac{\Delta^{4}}{2^{4}\cdot
24}&=\frac{\Delta^{2}}{8}\left(1-\frac{\Delta^{2}}{2\cdot 24}\right)\\\
&\geq\frac{\Delta^{2}}{8}\left(1-\frac{\pi^{2}}{2\cdot 24}\right).\end{split}$
(20)

Consequently, $\Delta=\Theta(\delta^{1/2})$, and this is the reason why
Quantum Walks are quadratically faster than their classical counterparts.

##  Appendix B  Mathematical description of the out-of-equilibrium quantum
Metropolis algorithm

In the appendix A we have reviewed the Szegedy quantum walk. In this appendix,
we present a quantum Metropolis-Hasting algorithm based on the use of Szegedy
walks.

The objective of the Metropolis-Hastings algorithm is sampling from the Gibbs
distribution
$\rho^{\beta}=\ket{\pi^{\beta}}\bra{\pi^{\beta}}=Z^{-1}(\beta)\sum_{\bm{\phi}\in\Omega}e^{-\beta
E(\bm{\phi})}\ket{\bm{\phi}}\bra{\bm{\phi}}$, where $E(\bm{\phi})$ is the
energy of a given configuration of angles of the molecule, $\beta$ plays the
role of the inverse of a temperature that will be lowered during the process,
and $Z(\beta)=\sum_{\bm{\phi}\in\Omega}e^{-\beta E(\bm{\phi})}$ a
normalization factor. $\Omega$ represents the configuration space, in our case
the possible values the torsion angles may take. One can immediately notice
that if $\beta$ is sufficiently large only the configurations with the lowest
possible energy will appear in when sampling from that state, with high
probability. Thus, we wish to prepare one such state to be able find the
configuration with the lowest energy, in our case the folded state of the
protein in nature.

One way to construct such distribution $\pi^{\beta}$ is to create a rapidly
mixing Markov Chain that has as equilibrium distribution $\pi^{\beta}$. Such
Markov Chain is characterized, at a given $\beta$, by a transition matrix
$\mathcal{W}$ that induces a random walk over the possible states. That is,
$\mathcal{W}$ maps a given distribution $p$ to another
$p^{\prime}=\mathcal{W}p$. Let us introduce some useful definitions: an
aperiodic walk is called irreducible if any state in $\Omega$ can be accessed
from any other state in $\Omega$, although not necessarily in a single step.
Additionally, we will say that a walk is reversible when it fulfills the
detailed balance condition

$\mathcal{W}_{j,i}^{\beta}\pi^{\beta}_{i}=\mathcal{W}_{i,j}^{\beta}\pi^{\beta}_{j}.$
(21)

Metropolis-Hastings algorithm uses the following transition matrix

$\mathcal{W}_{ij}=\begin{cases}T_{ij}A_{ij},&\text{if $i\neq j$}\\\
1-\sum_{k\neq j}T_{kj}A_{kj},&\text{if $i=j$},\end{cases}$ (22)

where, as given in the main text equation (1),

$A_{ij}=\min\left(1,e^{-\beta(E_{i}-E_{j})}\right),$ (23)

and

$T_{ij}=\begin{cases}\frac{1}{N}&\text{there is a move connecting $j$ to
$i$}\\\ 0,&\text{else},\end{cases}$ (24)

for $N$ the number of possible outgoing movements from state $j$, which we
assume to be independent of the current state, as is the case for our
particular problem. In the case of the Metropolis-Hastings algorithm the
detailed balance condition is fulfilled with the above definitions.

Having defined the Metropolis algorithm, we now want to quantize it using the
previous section. There have been several proposals to quantize the Metropolis
algorithm, as can be seen in the Appendix B. They often rely on slow evolution
from $\ket{\pi_{t}}\rightarrow\ket{\pi_{t+1}}$ using variations of Amplitude
Amplification, until a large $\beta$ has been achieved.

However, most often this Metropolis algorithms are used outside of
equilibrium, something that is not often worked out in the previous
approaches. There are at least two ways to perform the quantization of the
out-of-equilibrium Metropolis-Hastings algorithm lemieux2019efficient . The
first one, called ‘Zeno with rewind’ temme2011quantum  uses a simpler version
of previous quantum walks, where instead of amplitude amplification, quantum
phase estimation is used on operators $W_{j}$, so that measuring the first
eigenvalue means that we have prepared the corresponding eigenvector
$\ket{\psi_{t}^{0}}=\ket{\pi_{t}}$ somma2007quantum . Since the eigenvalue gap
is $\Delta_{t}$, the cost of Quantum Phase Estimation is $O(\Delta_{t}^{-1})$.

Define measurements $Q_{t}=\ket{\pi_{t}}\bra{\pi_{t}}$ and
$Q^{\perp}_{t}=\mathbf{1}-\ket{\pi_{t}}\bra{\pi_{t}}$, for $t$ an index
indicating the step of the cooling schedule. Performing these measurements
consists, as we have mentioned, in performing phase estimation of the
corresponding operator $\mathcal{W}_{t}$, at cost $\Delta_{t}^{-1}$. If a
measurement of type $Q^{\perp}_{t}$ indicates a restart of the whole algorithm
it is called ‘without rewind’. However, if measurement $Q^{\perp}_{t}$ is
obtained, one can perform phase estimation of $\mathcal{W}_{t-1}$. If
$|\braket{\pi_{t}}{\pi_{t-1}}|^{2}=F_{t}^{2}$, then the transition probability
between $Q_{j}\leftrightarrow Q_{t-1}$ and between
$Q^{\perp}_{t}\leftrightarrow Q^{\perp}_{t-1}$ is given by $F_{t}^{2}$; and
the transition probability between $Q_{t}\leftrightarrow Q^{\perp}_{t-1}$ and
between $Q^{\perp}_{t}\leftrightarrow Q_{t-1}$ is given by $1-F_{t}^{2}$, so
in a logarithmic number of steps one can recover state $\ket{\pi_{t-1}}$ or
$\ket{\pi_{t}}$.

The second proposal is to perform the unitary heuristic procedure
lemieux2019efficient

$\ket{\psi(L)}=W_{L}...W_{1}\ket{\pi_{0}}.$ (25)

This is in some ways the simplest way one would think of quantizing the
Metropolis-Hastings algorithm, implementing a quantum walk instead of a random
walk. It is very similar to the classical way of performing many steps of the
classical walk, slowly increasing $\beta$ until one arrives to the aim
temperature. In addition, this procedure is significantly simpler than the
previously explained ones because it does not require to perform phase
estimation on $W_{t}$.

Two more innovations are introduced by lemieux2019efficient . In the first
place, a heuristic Total Time to Solution (TTS) is defined. Assuming some
start distribution if operators $W$ are applied $t$ times, the probability of
success is given by $p(t)$. In order to be successful with constant
probability $1-\delta$ that means repeating the procedure
$\log(1-\delta)/\log(1-p(t))$ times. The total expected time to success is
then

$TTS(t):=t\frac{\log(1-\delta)}{\log(1-p(t))},$ (26)

as also indicated in the main text equation (7). For the unitary procedure

$TTS(L)=L\frac{\log(1-\delta)}{\log(1-|\braket{\pi^{g}}{W_{L}...W_{1}}{\pi_{0}}|^{2})}$
(27)

with $\pi^{g}$ the ground state of the Hamiltonian.

Additionally lemieux2019efficient  constructed an alternative to Szegedy
operator $\tilde{W}$ using a Boltzmann coin, such that the quantum walk
operator operates on three registers $\ket{x}_{S}\ket{z}_{M}\ket{b}_{C}$. $S$
stands for the codification of the state, $M$ for the codification of possible
movements, and $C$ the Boltzmann coin. Their operator $\tilde{W}$ is
equivalent to the Szegedy operator under a conjugation operator $Y$ that maps
moves to states and viceversa:

$\tilde{W}=RV^{\dagger}B^{\dagger}FBV,$ (28)

with,

$V:\ket{0}_{M}\rightarrow\sum_{j}N^{-1/2}\ket{j}$ (29a)
$\begin{split}B:&\ket{x}_{S}\ket{j}_{M}\ket{0}_{C}\rightarrow\\\
&\ket{x}_{S}\ket{j}_{M}\left(\sqrt{1-A_{x\cdot
z_{j},x}}\ket{0}+\sqrt{A_{x\cdot z_{j},x}}\ket{1}\right)\end{split}$ (29b)
$F:\ket{x}_{S}\ket{j}_{M}\ket{b}_{C}\rightarrow\ket{x\cdot
z^{b}_{j}}_{S}\ket{j}_{M}\ket{b}_{C}$ (29c)
$\begin{split}R:&\ket{0}_{M}\ket{0}_{C}\rightarrow-\ket{0}_{M}\ket{0}_{C}\\\
&\ket{j}_{M}\ket{b}_{C}\rightarrow\ket{j}_{M}\ket{b}_{C},\quad(j,b)\neq(0,0)\end{split}$
(29d)

Here $V$ proposes different moves, $B$ prepares the Boltzmann coin, $F$ flips
the bits necessary to prepare the new state, conditional on the Boltzmann coin
being in state 1, and $R$ is a reflection operator on state $(0,0)$ for the
coin and movement registers. Although our encoding of the operators and states
is slightly different, the algorithm that we have used is this one, mainly due
to its simplicity.

## References

  * (1) U. Consortium, “Uniprot: a worldwide hub of protein knowledge,” Nucleic Acids Research, vol. 47, no. D1, pp. D506–D515, 2019.
  * (2) H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov, and P. E. Bourne, “The protein data bank,” Nucleic Acids Research, vol. 28, no. 1, pp. 235–242, 2000.
  * (3) N. Perdigão, J. Heinrich, C. Stolte, K. S. Sabir, M. J. Buckley, B. Tabor, B. Signal, B. S. Gloss, C. J. Hammang, B. Rost, et al., “Unexpected features of the dark proteome,” Proceedings of the National Academy of Sciences, vol. 112, no. 52, pp. 15898–15903, 2015.
  * (4) A. Bhowmick, D. H. Brookes, S. R. Yost, H. J. Dyson, J. D. Forman-Kay, D. Gunter, M. Head-Gordon, G. L. Hura, V. S. Pande, D. E. Wemmer, et al., “Finding our way in the dark proteome,” Journal of the American Chemical Society, vol. 138, no. 31, pp. 9730–9742, 2016.
  * (5) N. Perdigão and A. Rosa, “Dark proteome database: studies on dark proteins,” High-throughput, vol. 8, no. 2, p. 8, 2019.
  * (6) P. N. Bryan and J. Orban, “Proteins that switch folds,” Current Opinion in Structural Biology, vol. 20, no. 4, pp. 482–488, 2010.
  * (7) A. K. Dunker, J. D. Lawson, C. J. Brown, R. M. Williams, P. Romero, J. S. Oh, C. J. Oldfield, A. M. Campen, C. M. Ratliff, K. W. Hipps, et al., “Intrinsically disordered protein,” Journal of Molecular Graphics and Modelling, vol. 19, no. 1, pp. 26–59, 2001.
  * (8) R. Das and D. Baker, “Macromolecular modeling with rosetta,” Annu. Rev. Biochem., vol. 77, pp. 363–382, 2008.
  * (9) University of Washington<EMAIL_ADDRESS>boinc.bakerlab.org, 2021.
  * (10) R. Das, B. Qian, S. Raman, R. Vernon, J. Thompson, P. Bradley, S. Khare, M. D. Tyka, D. Bhat, D. Chivian, et al., “Structure prediction for casp7 targets using extensive all-atom refinement with rosetta@ home,” Proteins: Structure, Function, and Bioinformatics, vol. 69, no. S8, pp. 118–128, 2007.
  * (11) W. E. Hart and S. Istrail, “Robust proofs of np-hardness for protein folding: general lattices and energy potentials,” Journal of Computational Biology, vol. 4, no. 1, pp. 1–22, 1997.
  * (12) B. Berger and T. Leighton, “Protein folding in the hydrophobic-hydrophilic (hp) is np-complete,” in Proceedings of the Second Annual International Conference on Computational Molecular Biology, pp. 30–39, 1998.
  * (13) A. Kryshtafovych, T. Schwede, M. Topf, K. Fidelis, and J. Moult, “Critical assessment of methods of protein structure prediction (casp)—round xiii,” Proteins: Structure, Function, and Bioinformatics, vol. 87, no. 12, pp. 1011–1020, 2019.
  * (14) A. Senior, R. Evans, J. Jumper, J. Kirkpatrick, L. Sifre, T. Green, C. Qin, A. Zidek, A. Nelson, A. Bridgland, et al., “Improved protein structure prediction using potentials from deep learning,” Nature, 2020.
  * (15) M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek, and G. R. Hutchison, “Avogadro: an advanced semantic chemical editor, visualization, and analysis platform,” Journal of Cheminformatics, vol. 4, no. 1, p. 17, 2012.
  * (16) P. Wocjan and A. Abeyesinghe, “Speedup via quantum sampling,” Physical Review A, vol. 78, no. 4, p. 042336, 2008.
  * (17) R. Somma, S. Boixo, and H. Barnum, “Quantum simulated annealing,” arXiv preprint arXiv:0712.1008, 2007.
  * (18) R. D. Somma, S. Boixo, H. Barnum, and E. Knill, “Quantum simulations of classical annealing processes,” Physical Review Letters, vol. 101, no. 13, p. 130504, 2008.
  * (19) K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete, “Quantum metropolis sampling,” Nature, vol. 471, no. 7336, p. 87, 2011\.
  * (20) M.-H. Yung and A. Aspuru-Guzik, “A quantum–quantum metropolis algorithm,” Proceedings of the National Academy of Sciences, vol. 109, no. 3, pp. 754–759, 2012.
  * (21) J. Lemieux, B. Heim, D. Poulin, K. Svore, and M. Troyer, “Efficient Quantum Walk Circuits for Metropolis-Hastings Algorithm,” Quantum, vol. 4, p. 287, June 2020.
  * (22) M. Szegedy, “Quantum speed-up of markov chain based algorithms,” in 45th Annual IEEE Symposium on Foundations of Computer Science, pp. 32–41, IEEE, 2004\.
  * (23) R. Babbush, A. Perdomo-Ortiz, B. O’Gorman, W. Macready, and A. Aspuru-Guzik, “Construction of energy functions for lattice heteropolymer models: a case study in constraint satisfaction programming and adiabatic quantum optimization,” arXiv preprint arXiv:1211.3422, 2012.
  * (24) A. Robert, P. K. Barkoutsos, S. Woerner, and I. Tavernelli, “Resource-efficient quantum algorithm for protein folding,” arXiv preprint arXiv:1908.02163, 2019.
  * (25) A. Perdomo-Ortiz, N. Dickson, M. Drew-Brook, G. Rose, and A. Aspuru-Guzik, “Finding low-energy conformations of lattice protein models by quantum annealing,” Scientific Reports, vol. 2, p. 571, 2012.
  * (26) M. Fingerhuth, T. Babej, et al., “A quantum alternating operator ansatz with hard and soft constraints for lattice protein folding,” arXiv preprint arXiv:1810.13411, 2018.
  * (27) T. Babej, M. Fingerhuth, et al., “Coarse-grained lattice protein folding on a quantum annealer,” arXiv preprint arXiv:1811.00713, 2018.
  * (28) A. Perdomo, C. Truncik, I. Tubert-Brohman, G. Rose, and A. Aspuru-Guzik, “Construction of model hamiltonians for adiabatic quantum computation and its application to finding low-energy conformations of lattice protein models,” Physical Review A, vol. 78, no. 1, p. 012320, 2008.
  * (29) C. Outeiral, G. M. Morris, J. Shi, M. Strahm, S. C. Benjamin, and C. M. Deane, “Investigating the potential for a limited quantum speedup on protein lattice problems,” arXiv preprint arXiv:2004.01118, 2020.
  * (30) V. K. Mulligan, H. Melo, H. I. Merritt, S. Slocum, B. D. Weitzner, A. M. Watkins, P. D. Renfrew, C. Pelissier, P. S. Arora, and R. Bonneau, “Designing peptides on a quantum computer,” bioRxiv, p. 752485, 2020.
  * (31) L. Banchi, M. Fingerhuth, T. Babej, C. Ing, and J. M. Arrazola, “Molecular docking with gaussian boson sampling,” Science Advances, vol. 6, no. 23, p. eaax1950, 2020.
  * (32) The $\alpha$-helix and the $\beta$-sheet correspond to two common structures found in protein folding. Such structures constitute what is called the secondary structure of the protein, and are characterised because $(\phi,\psi)=(-\pi/3,-\pi/4)$ in the $\alpha$-helix, and $(\phi,\psi)=(-3\pi/4,-3\pi/4)$ in the $\beta$-sheet, due to the hydrogen bonds that happen between backbone amino groups NH and backbone carboxy groups CO.
  * (33) R. Von Mises, Mathematical Theory of Probability and Statistics.  Academic Press, 2014.
  * (34) J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. A. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke, M. L. Abrams, et al., “Psi4: an open-source ab initio electronic structure program,” Wiley Interdisciplinary Reviews: Computational Molecular Science, vol. 2, no. 4, pp. 556–565, 2012.
  * (35) L. K. Grover, “Quantum mechanics helps in searching for a needle in a haystack,” Physical Review Letters, vol. 79, no. 2, p. 325, 1997.
  * (36) F. Magniez, A. Nayak, J. Roland, and M. Santha, “Search via quantum walk,” SIAM Journal on Computing, vol. 40, no. 1, pp. 142–164, 2011.
  * (37) T. Albash and D. A. Lidar, “Demonstration of a scaling advantage for a quantum annealer over simulated annealing,” Physical Review X, vol. 8, no. 3, p. 031016, 2018.
  * (38) E. Alcaide, “Minifold: a deeplearning-based mini protein folding engine.” https://github.com/EricAlcaide/MiniFold/, 2019.
  * (39) H. Abraham et al., “Qiskit: An open-source framework for quantum computing,” 2019.
  * (40) Amazon.com, Inc., “Amazon web services.” aws.amazon.com, 2021.
  * (41) S. Kim, J. Chen, T. Cheng, A. Gindulyte, J. He, S. He, Q. Li, B. A. Shoemaker, P. A. Thiessen, B. Yu, et al., “Pubchem 2019 update: improved access to chemical data,” Nucleic Acids Research, vol. 47, no. D1, pp. D1102–D1109, 2019.
  * (42) F. Jensen, “Atomic orbital basis sets,” Wiley Interdisciplinary Reviews: Computational Molecular Science, vol. 3, no. 3, pp. 273–295, 2013.
  * (43) T. Helgaker, P. Jorgensen, and J. Olsen, Molecular electronic-structure theory.  John Wiley & Sons, 2014.
  * (44) L. Grover and T. Rudolph, “Creating superpositions that correspond to efficiently integrable probability distributions,” arXiv preprint quant-ph/0208112, 2002.
  * (45) S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 671–680, 1983.
  * (46) K. Temme, S. Bravyi, and J. M. Gambetta, “Error mitigation for short-depth quantum circuits,” Physical Review Letters, vol. 119, no. 18, p. 180509, 2017.
  * (47) R. LaRose, A. Mari, P. J. Karalekas, N. Shammah, and W. J. Zeng, “Mitiq: A software package for error mitigation on noisy quantum computers,” 2020.
  * (48) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al., “Tensorflow: A system for large-scale machine learning,” in 12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16), pp. 265–283, 2016\.
  * (49) F. Chollet et al., “Keras.” https://github.com/fchollet/keras, 2015\.
  * (50) R. S. Sutton and A. G. Barto, Reinforcement learning: An introduction.  MIT press, 2018.
  * (51) P. J. Van Laarhoven and E. H. Aarts, “Simulated annealing,” in Simulated annealing: Theory and applications, pp. 7–15, Springer, 1987.
  * (52) Student, “The probable error of a mean,” Biometrika, pp. 1–25, 1908.
  * (53) R. Portugal, Quantum walks and search algorithms.  Springer, 2013.
  * (54) A. Ambainis, “Quantum walk algorithm for element distinctness,” SIAM Journal on Computing, vol. 37, no. 1, pp. 210–239, 2007.
  * (55) G. D. Paparo and M. Martin-Delgado, “Google in a quantum network,” Scientific Reports, vol. 2, p. 444, 2012.
  * (56) G. D. Paparo, M. Müller, F. Comellas, and M. A. Martin-Delgado, “Quantum google in a complex network,” Scientific Reports, vol. 3, p. 2773, 2013\.
  * (57) G. D. Paparo, V. Dunjko, A. Makmal, M. A. Martin-Delgado, and H. J. Briegel, “Quantum speedup for active learning agents,” Physical Review X, vol. 4, no. 3, p. 031002, 2014.
  * (58) Notice that in most texts the definition of $M$ does not explicitly include $S$. It is assumed implicitly though.