# An Efficient MCMC Approach to Energy Function Optimization
in Protein Structure Prediction

Lakshmi A. Ghantasala<EMAIL_ADDRESS>Corresponding Author Elmore Family
School of Electrical and Computer Engineering, Purdue University  Risi Jaiswal
Elmore Family School of Electrical and Computer Engineering, Purdue University
Supriyo Datta  Elmore Family School of Electrical and Computer Engineering,
Purdue University

###### Abstract

Protein structure prediction is a critical problem linked to drug design,
mutation detection, and protein synthesis, among other applications. To this
end, evolutionary data has been used to build contact maps which are
traditionally minimized as energy functions via gradient descent based schemes
like the L-BFGS algorithm. In this paper we present what we call the
Alternating Metropolis-Hastings (AMH) algorithm, which (a) significantly
improves the performance of traditional MCMC methods, (b) is inherently
parallelizable allowing significant hardware acceleration using GPU, and (c)
can be integrated with the L-BFGS algorithm to improve its performance. The
algorithm shows an improvement in the energy of 8.17% to 61.04% (average
38.9%) over traditional MH and 0.53% to 17.75% (average 8.9%) over traditional
MH with intermittent noisy restarts, tested across 9 proteins from recent CASP
competitions. We go on to map the Alternating MH algorithm to a GPGPU which
improves sampling rate by 277x and improves simulation time to a low energy
protein prediction by 7.5x to 26.5x over CPU. We show that our approach can be
incorporated into state-of-the-art protein prediction pipelines by applying it
to trRosetta2’s energy function and also the distogram component of
Alphafold1’s energy function. Finally, we note that specially designed
probabilistic computers (or p-computers) can provide even better performance
than GPU’s for MCMC algorithms like the one discussed here.

Markov chain monte carlo, Metropolis-Hastings, protein structure prediction,
optimization

††preprint: APS/123-QED

##  I Introduction

The problem of protein folding is a decades long quest launched in large part
by Anfinsen’s thermodynamic hypothesis that a protein’s geometry can be
inferred from only its’ amino acid sequence [1]. In the last decade, co-
evolutionary sequence data for a given protein has been used to build contact
maps from which the 3-D geometry of a protein can be inferred. Contact maps
constitute a set of probabilities for inter-residue distances or orientations.
The structures that satisfy these probabilities most are then the predictions
for the 3-D geometry of the protein for which the contact map is intended [2,
3]. While evolutionary couplings may be statistically derived [4, 5], deep
neural networks have seen resounding success in predicting a protein’s contact
map [6, 7, 8]. Deep neural networks were notably used by the Alphafold1 team
to come in first in CASP13 (2018) [9], and again by trRosetta2 in CASP14
(2020).

Predicting highly accurate contact maps from co-evolutionary data, i.e
optimizing inter-residue distance and orientation probabilities, is the first
half of the problem. These probabilities may then be formulated into an energy
function whose minimum describes the optimal conformation of the protein. The
goal is to find a set of residue ($C_{\beta}$) positions that minimize the
derived energy function and related force-fields. For both Alphafold1 and
trRosetta2, non-monotone L-BFGS [10], a second order gradient descent method,
was used. On the other hand, it has recently been shown that sampling
approaches may be better suited to find optimal states in non-convex energy
landscapes when compared to traditional gradient descent approaches [11]. This
naturally raises the suspicion that a direct sampling approach to minimizing
contact map-based energy functions may lead to better results than gradient
descent.

We present a novel sampling approach termed Alternating Metropolis-Hastings
(Alternating MH) that acts directly on the torsion angles (ex. $\varphi$,
$\psi$) of a protein to find optimal conformations. Alternating MH improves
upon traditional MH by introducing local energies by which proposals can be
judged. Localized MH limits the scope of the energy to short segments within
the full protein. By limiting the scope of the energy function, we essentially
assign a ‘local’ energy to a proposal state that is based only on the
immediate surroundings of the modified variable. The current energy is also
determined on the same range, and the comparison now judges whether the
proposal improves local energy or not. Alternating MH ’alternates’ between
local and global MH to optimize an energy function.

It should be noted that recent works like Alphafold 2 [12] and RoseTTAFold
[13] based on end-to-end prediction of protein structure using transformer
networks have observed remarkable success. These machine learning pipelines
have made significant improvements over the previous contact-map based methods
that our Alternating MH is shown to expedite. However, our method should be
adaptable to new state-of-the-art techniques. Indeed it may be even more
useful in future for $\it{in\ vivo}$ protein structure prediction which will
require sampling from stochastic environments as opposed to the controlled
environments encountered $\it{in\ vitro}$.

##  II Methods

###  II.1 Alternating Metropolis-Hastings Algorithm

Traditional Metropolis-Hastings (MH) forms the first stage of Alternating MH.
Proposal states are repeatedly generated, beginning at an input initial state,
and sequentially accepted or rejected according to an acceptance probability
($P_{accept}$). The proposal scheme acts on an $n$-variable state at step $i$,
$s^{i}_{1,\ldots,n}$ , where $s^{i}_{even}=\varphi^{i}_{1,\ldots,r}$ and
$s^{i}_{odd}=\psi^{i}_{1,\ldots,r}$ for a protein with $r$ residues. This
mapping, used by Alphafold1, is illustrated in Fig. 2. While different
implementations may use different mappings (ex. trRosetta2 follows a 5 angle
parameterization), they may all be mapped to the $s$ vector, which is finally
used to generate proposals. The proposal mechanism is a very simple ’single-
flip’ mechanism where a Gaussian noise is added to a single torsion angle in
$s$ per proposal. After $n$ proposals, we arrive at an entirely new sample.
Energy is calculated as a function of $s$ and the coordinates of atoms in the
resulting structures. The probability of acceptance of a proposal for
traditional MH is described in Eq. 1a:

$\displaystyle P_{accept}$ $\displaystyle(s^{i}_{1,\ldots,n})=$ (1a)
$\displaystyle\min\Big{[}1,\frac{P(s^{i+1}_{1,\ldots,n})P(s^{i+1}_{1,\ldots,n}|s^{i}_{1,\ldots,n})}{P(s^{i}_{1,\ldots,n})P(s^{i}_{1,\ldots,n}|s^{i+1}_{1,\ldots,n})}\Big{]}$

where superscripts define the iteration of the state and subscripts define the
range of variables being considered within the state. The subscripts will
become relevant when a localized version of this algorithm is defined. We can
simplify the acceptance criterion to remove dependence on the conditional
probability, as a symmetric, Gaussian transition distribution is used:

$\displaystyle
P(s^{i+1}_{1,\ldots,n}|s^{i}_{1,\ldots,n})=N(s^{i}_{k},\sigma).$ (1b)

The marginal component of the acceptance criterion may be defined as:

$\displaystyle P(s^{i}_{1,\ldots,n})=e^{-V(x)}.$ (1c)

where $x=G(s^{i}_{1,\ldots,n})$ defines the coordinates ($x$) for all
$C_{\beta}$ atoms in the protein as a function of the torsion angles of the
chain.

Alternating MH is a combination of two stages. Stage one is traditional MH as
described in Eqs. 1. Stage two is termed Localized Metropolis-Hastings (see
Eqs. 2), which uses a modified energy function that assigns energy to a
proposal based only on some local segment of $s$. The acceptance of a proposal
is now based on the energy of the segment of the residue chain:

$\displaystyle
P_{accept}(s^{i}_{1,\ldots,n})=\min\Big{[}1,\frac{P(s^{i+1}_{k-d,\ldots,k+d})}{P(s^{i}_{k-d,\ldots,k+d})}\Big{]}$
(2a)

where $k$ is the index of the variable being modified within the state and $d$
is the radius of the segment. To get the marginal probability of a segment of
a proposal, a new localized potential must be defined which returns the
segment’s energy:

$\displaystyle P(s^{i}_{k-d,\ldots,k+d})=e^{-V_{local}(x,k,d)}$ (2b)

Generally, $V$ requires summing over all inter-residue (and inter-atomic)
energies for a protein; we may define a localized energy function
$V_{\textrm{local}}$ by summing over a subset of the inter-residue energies
defined by $k$ and $d$.

Figure 1:  The Alternating Mh scheme cycles between generating $N_{global}$
samples via traditional MH, and $N_{local}$ samples via a novel localized MH.
The traditional MH may be thought of as the global optimization component of
the process, akin to L-BFGS, while the localized MH may be compared to ’noisy
restarts’, an intelligent way to escape from local minima.

Localized MH follows the same Gaussian transition probability defined in Eq.
1b. Alternating MH alternates between local and global MH, judging $N_{local}$
proposals according to local energies, then $N_{global}$ proposals according
to global energy, as in Fig. 1. The pseudo code for the Alternating MH scheme
is shown in Alg. 1.

input : $i=0,d=10,s^{0}=\\{180^{\circ},\ldots,180^{\circ}\\}$

output : $s\leftarrow\textrm{optimal angles}$

for  _$l\leftarrow 0$ to $cycles$_ do

/* Traditional Metropolis-Hastings */

for  _$i\leftarrow 1$ to $N_{\textrm{global}}$_ do

for  _$k\leftarrow 1$ to length($s$)_ do

$s^{{}^{\prime}}\leftarrow s^{i-1}$

$s^{{}^{\prime}}_{k}\leftarrow s^{i-1}_{k}+N(0,\sigma)$

$P_{\textrm{accept}}=\min\Big{[}1,\frac{e^{-V(\textrm{G}(s^{i-1}))}}{e^{-V(\textrm{G}(s^{{}^{\prime}}))}}\Big{]}$

if  _$\textrm{rand}(0,1) <P_{\textrm{accept}}$_ then

$s^{i}\leftarrow s^{{}^{\prime}}$

else

$s^{i}\leftarrow s^{i-1}$

end if

end for

end for

/* Localized Metropolis-Hastings */

for  _$i\leftarrow 1$ to $N_{\textrm{local}}$_ do

for  _$k\leftarrow 1$ to length($s$)_ do

$s^{{}^{\prime}}\leftarrow s^{i-1}$

$s^{{}^{\prime}}_{k}\leftarrow s^{i-1}_{k}+N(0,\sigma)$

$P_{\textrm{accept}}=\min\Big{[}1,\frac{e^{-V_{\textrm{local}}(\textrm{G}(s^{i-1}),k,d)}}{e^{-V_{\textrm{local}}(\textrm{G}(s^{{}^{\prime}}),k,d)}}\Big{]}$

if  _$\textrm{rand}(0,1) <P_{\textrm{accept}}$_ then

$s^{i}\leftarrow s^{{}^{\prime}}$

else

$s^{i}\leftarrow s^{i-1}$

end if

end for

end for

end for

Algorithm 1 Alternating Metropolis-Hastings Figure 2:  (A) The mapping of an
unfolded protein (T0955 from CASP13) to variables is shown here. The dihedral
angles $(\varphi,\psi)$ are captured by a single vector s such that
$s^{i}_{\textrm{even}}=\varphi^{i}_{1,\ldots,r}$ and
$s^{i}_{\textrm{odd}}=\psi^{i}_{1,\ldots,r}$. The uncolored atoms attached to
the $C_{\beta}$ (blue) atoms are the side chains for each residue. The
unfolded protein in (A) is the initial ‘current state’ in the following flow
chart. (B) depicts the flow chart for an MH based minimization process.
Alternating MH impacts step 3, where it modifies the energy used for the MH
criterion between traditional MH and localized MH.

####  II.1.1 Localization Radius Parameter $d$

The localization radius parameter describes the number of variables in $s$
from variable $k$ that we may consider for the MH criterion. This gives us a
powerful knob with which to tune how the MCMC procedure folds the protein. The
smaller $d$ is, the fewer residues each proposed fold has to satisfy, and the
larger the acceptance rate is. Traditional MH, by this metric, has the
smallest acceptance ratio where every proposal must satisfy, to some degree,
every residue in the chain.

As the protein folding problem may be captured by a 1-D state array ($s$), $d$
amounts to looking to the left and right of the $k^{th}$ variable in $s$. If
the state vector is 2 (or more)-Dimensional, the segment may grow as a circle
($n$-d sphere), radiating outwards from the center $k$ variable with radius
$d$. The advantage then of defining a radius parameter, rather than a range
parameter, is in easily extrapolating this algorithm to problems with higher
dimensional state vectors.

One can cascade multiple localization ranges to consider folds of multiple
sizes. While results for this work alternate between global MH, i.e $d=2n$
where $n$ is the length of $s$, and a localized MH with $d=10$, further ranges
could also be included. Intuitively, it would seem the localization ranges
should reflect the sizes of local structures in the problem, though future
work will have to determine how best to cater the radius to the problem at
hand.

###  II.2  Applying Alternating MH to Protein Structure Prediction

Two implementations of the Alternating MH scheme are demonstrated for contact
map-based energy minimization. The first implementation (described by Eq. 3)
minimizes an energy function built with Alphafold1’s published distograms [2,
14]. The energy function is as follows:

$\displaystyle V(x)=-\sum_{i,j,i\neq j}\log P(|x_{i}-x_{j}||S,MSA(S))$ (3a)

Where $x=G(\varphi,\psi)$, $S$ is the sequence of amino acids to be folded,
and $MSA(S)$ is the Multiple Sequence Alignment of that sequence. This
function is a part of the total energy used in the Alphafold1 pipeline, which
also includes a torsion potential and a physics based term that includes van-
der-waals forces; these additional terms are discounted in this paper to make
an accelerated GPGPU implementation more reliable. To find the energy for a
localized segment, the log probabilities are summed over only within the
desired range:

$\displaystyle V_{\textrm{local}}(x,k,d)=$
$\displaystyle-\sum^{\lfloor\frac{k}{2}\rfloor+\lfloor\frac{d}{2}\rfloor}_{i,j=\lfloor\frac{k}{2}\rfloor-\lfloor\frac{d}{2}\rfloor,i\neq
j}\log P(|x_{i}-x_{j}||S,MSA(S))$ (3b)

where $k$ is the index of the variable modified between subsequent states, and
$d$ is the radius of the segment within the angle vector ($s$). As there are 2
torsion angles per residue in the parameterization described by Alphafold1,
the angle being modified is mapped to its corresponding $C_{\beta}$ atom in
the chain by dividing by 2 and flooring.

Details regarding the GPGPU implementation of the minimization of this
distogram-only energy function are in section II.3.

The second implementation, described in Eq. II.2 minimizes an energy function
that is a part of the trRosetta2 pipeline [3]. Unlike with Alphafold1, the
complete energy function was minimized with Alternating MH, which includes
both restraint and physical forces:

$\displaystyle V(s)=$ $\displaystyle
V_{\textrm{distance}}(G(s))+V_{\textrm{torsion}}(s)+$ $\displaystyle
V_{\textrm{rosetta}}(G(s))$

where $V_{\textrm{rosetta}}$ includes energies corresponding to ramachandran
(rama), hydrogen bonding (cen_hb), and van der Waals (vdw) forces [3]. A
localized version of this function was derived via the $get\\_sub\\_score$
function in pyrosetta. The trRosetta2 pipeline parameterizes $C_{\beta}$ atoms
according to $(\varphi,\theta,\omega)$, which is different than Alphafold’s
parameterization, though once the information is captured in the angle vector
($s$), the remaining procedure is unchanged. Proposals are still generated and
judged via $s$.

The original minimization process for trRosetta2 involves modifying the
weights of the score function as the minimization scheme iterates; for the
sake of an exact comparison between L-BFGS and Alternating MH, the process was
simplified to minimizing a single, constantly weighted energy function.
Furthermore, the final relaxation procedure (FastRelax) is not considered. The
final RMSD values, then, will not match with trRosetta2’s published results,
but should make for a robust comparison with the Alternating MH scheme.

In both implementations, $s$ is initialized to {180°}, as shown in Fig. 2. The
proposal scheme picks $s^{i}_{k}$ to perturb by a 0 mean Gaussian noise at
each proposal, incrementing $k$ from 0 to $n$ where $n$ is the length of $s$.
The variance of the Gaussian ($\sigma$) is gradually reduced over the course
of an Alternating MH run; for this work, $\sigma=10\rightarrow 4$. This
ensures that initial proposals result from larger swings in $s$, and later
proposals become tighter. The full Alt-MH procedure was repeated 10 times per
tested protein, where each restart would begin at the lowest energy found
during the previous run. If a lower energy was not found during the previous
Alt-MH run, the range of $\sigma$ changes to $10/L\rightarrow 1/L$ where $L$
is the number of runs where a lower energy was not found.

Figure 3:  (A) The Alternating Metropolis-Hastings method presented in this
paper (column 3) is benchmarked against various optimization schemes that
include MCMC, L-BFGS, and mixed schemes. Each scheme is applied 10 times to
each of 9 different proteins sampled from the previous 3 CASP competitions;
the minimum energies observed over each run are plotted. All runs are for
10000 samples. While the L-BFGS implementation in column 4 is based
specifically on trRosetta2’s implementation, Alphafold1 and trRosetta1 also
use a variation of this approach, hence their inclusion in the label. (B) The
lowest energy prediction for T0974s1 by the standard L-BFGS scheme (blue) and
our proposed Alternating MH scheme (blue) shows Alt-MH is able to find a lower
energy state with improved RMSD to the original structure (green).

###  II.3  Hardware Acceleration of Alternating MH

The Metropolis Hastings algorithm is generally amenable to hardware that
allows multiple chains to be run in parallel. The more chains the hardware
accommodates, the larger the space that is searched and consequently the
better the chances of finding lower energy ground states.

A General Purpose Graphics Processing Unit (GPGPU) was chosen as a suitable
hardware accelerator for Alt-MH, as it is a parallel computing platform that
naturally accommodates the complex floating-point operations required for this
problem (ex. cos, sin, log, sqrt), as opposed to an FPGA. As a note, we use
the terms GPGPU and GPU interchangeably in this manuscript as the Nvidia P100
may be seen as both. A simplified diagram of the GPU architecture as an
accelerator for the Alternating MH algorithm is presented in 4. The logical
blocks in the architecture of a p-computer [15], shown in 4.b, can be matched
1-to-1 with our mapping of AMH to a GPU.

Each block in the GPU houses an MCMC chain while an array of threads within
the block parallelize the operations within an MCMC chain. These two forms of
parallelism are key to achieving the presented performance.

####  II.3.1 Intra-Block parallelism

Alt-MH, and in-fact MCMC algorithms in general, involve(s) 4 steps; random
number generation, proposal generation, energy calculation, and criterion
evaluation; these steps may be parallelized to observe faster execution. For
the protein structure prediction problem, there are 3 steps that may take
advantage of intra-block parallelism:

1\. RNG: Random numbers are required for a) generating a proposal, and b)
evaluating the MH criterion. To generate these, each MCMC block in the GPU
maintains two LFSRs, assigned to the first 2 threads in each block. Each LFSR
produces a pseudo-random, uniformly distributed random number. One of these is
converted to a Gaussian-distributed random number via the Box-Muller method,
as required for proposal generation.

2\. Proposal Generation: A proposal protein structure is generated by adding a
Gaussian distributed, 0 mean value to each $\varphi$ or $\psi$ bond angle in
the protein. Every time a bond angle is modified, the subsequent (if the bond-
of-interest is phi) or preceding (if the bond-of-interest is psi) atoms must
rotate around the axis of the bond-of-interest. On CPU, this rotation is done
atom by atom; on GPU, each atom is assigned to its own thread, which allows
rotation of atoms in parallel. Atoms are rotated $n$ at a time, where $n$ is
the length of the angle vector $s$.

3\. Energy Calculation: The Alphafold1 implementation of Alt-MH relies on the
calculation of pairwise energies between residues; pairwise energies are then
summed across all ij residue pairs that satisfy
$\lfloor\frac{k}{2}\rfloor-\lfloor\frac{d}{2}\rfloor<i<j<\lfloor\frac{k}{2}\rfloor+\lfloor\frac{d}{2}\rfloor$,
where k is the residue housing the bond-of-interest and $d$ is the
localization range. In the GPU, each thread determines the pairwise energy of
a single ij pair. Each thread maintains a sum of its’ ij pairs until all
pairwise energies are accounted for. Then parallel reduction is employed to
efficiently sum the resultant energy terms across every thread in the block.

The effectiveness of parallelizing intra-block operations is shown in Fig 4.e.
As the number of threads grows, the simulation time of reaching 8.2e4
proposals reduces.

####  II.3.2 Inter-Block parallelism

Multiple Alt-MH blocks each with their own RNG will explore different slices
of the state space. The wider the search, the better the chances of finding
lower energy states. Accomodation of multiple parallel blocks also opens the
door for various coordinated MCMC schemes to enter the picture, ex. parallel
tempering [16], multiple-try metropolis [17], evolutionary monte-carlo [18],
and interchain adaptive MCMC [19] to name a few. These methods may work hand
in hand with Alt-MH. Alt-MH deals with improving a single MCMC block while
these methods work to coordinate multiple MCMC blocks to gain some advantage.

Figure 4:  (A) describes how the Alternating-MH algorithm may take advantage
of parallelized hardware to improve performance. The RNG and Kernel described
by the generic p-computer architecture in (B) can be found inside each MCMC
block. The operations within an MCMC block may be parallelized via threads in
the GPU as shown in (A). (C) measures the improvement in minimum energy found
via 1000 parallel blocks over 10 parallel blocks. While the mean (dotted
lines) of the two runs largely overlap, the min and max (solid lines) widen
between the 10 and 1000 block runs. (D) measures the effect of the number of
threads in a block on simulation time. (E) notes an overall improvement in
simulation time between the CPU (Intel I3-10100 @ 3.6Ghz) and GPU (Nvidia
P100) implementations across 9 protein domains, annotated with their CASP
labels.

##  III Results

The trRosetta2 based implementation is used as a benchmark to compare Alt-MH
against standard MCMC approaches, nonmonotone L-BFGS, and mixed approaches, as
seen in Fig. 3. Each optimization scheme was run for 10,000 samples. A sample
in the context of L-BFGS is when every variable in the state has been
perturbed according to its’ gradient once. A sample in the context of MH is
when every variable in the state has had a chance to be perturbed once. Note
that for MH based schemes, a sample need only be a full set of attempted
perturbations, as opposed to gradient descent based schemes where a sample is
a full set of perturbations.

The first column shows traditional Metropolis-Hastings with annealing, i.e
$0.1<\beta<1$ where $\beta$ is inverse temperature. The second column is
traditional MH, with the same annealing, but with a noisy restart every 100
samples. The noisy restarts add a uniform random number between -10 and 10 to
every torsion angle in the protein. Column 3 is the Alternating-MH algorithm
presented in this paper. The benchmark was done with 80 samples of global-MH
and 20 samples of localized MH with localization radius $d=5$. Perturbation
variance for proposals was reduced from $10^{\circ}$ to $4^{\circ}$ over the
course of each run, so later proposals are made within a shorter distance from
the current state than earlier proposals. This method improves upon the
energies of Traditional MH (column 1) by 8.17% for T0980s2 to 62.27% for
T1046s1 (and 61.04% for T1006) with average improvement at 38.9%. We also
observe improvement over MH with noisy restarts (column 2) of 0.53% for
T0953s1 to 17.75% for T1046s1.

The L-BFGS method in the benchmark is the lbfgs_armijo_nonmonotone minmover
from pyrosetta, with implementation details matching trRosetta2’s code.
Because this is a nonmonotone optimization scheme, it was benchmarked by
running in batches of 100 rather than in individual steps; performance drops
significantly when run one iteration at a time. Column 5, L-BFGS w/ localized
MH, uses the same L-BFGS function as column 4, but shortens the L-BFGS samples
to 80 and replaces the noisy restarts with 20 samples of localized MH, the
scheme described in this manuscript. Despite running for 20% fewer samples of
L-BFGS, column 5 exhibits predictions largely on par with the pure L-BFGS
method, i.e to within 3.5% for 7 of the 9 proteins, and within 12.8% worst
case. The last method, Metropolis-Hastings with L-BFGS generated proposals,
implements a Metropolis-adjusted L-BFGS algorithm in a similar fashion to the
well known Metropolis-adjusted Langevin algorithm (MALA) [20]. Proposals are
generated via lbfgs_armijo_nonmonotone run for 20 samples, and accepted or
rejected according to the metropolis-hastings criterion 10000 times. This
approach proved to be exactly on par with the standard L-BFGS approach (column
4), coming to within 1% worse than L-BFGS for 5 proteins (worst case 6.4% for
T1006), and beating L-BFGS by 1% for 2 proteins.

There is one odd-ball result that stands out amongst all the tested proteins,
T0974s1. Every scheme with some MH component finds a lower energy state for
T0974s1 compared to the standard L-BFGS minimization. Future work would have
to discover the quality that causes this. For T0974s1, our algorithm
outperforms the standard L-BFGS approach by 10.51% in energy and 32.5% in
RMSD, as seen in Fig. 3.B. L-BFGS with Localized MH (Column 4) outperforms the
standard approach by 23.73% in energy, and Metropolis-adjusted L-BFGS (Column
5) outperforms by 32.34%.

The Alternating-MH implementation described in Eq. 3 was mapped to an NVIDIA
P100 GPU running on Google Colaboratory, the results of which are shown in
Fig. 4. Fig.4.A describes how the algorithm fits into the architecture of a
GPU to take advantage of the hardware’s parallelism. The architecture of a
p-computer can be overlaid on the GPU architecture to show that this GPU
implementation fits into the Kernel-RNG paradigm for a p-computer described by
Kaiser et al [15, 21], shown in 4.B.

Fig 4.C showcases the benefits of having multiple parallel MCMC blocks; The
1000 parallel block run finds a structure with 12.97% lower energy than 10
blocks run in parallel; this directly translates into a predicted protein
structure with 23.0% lower RMSD to the original protein. Future works may
explore coordinated MCMC schemes that may lead to significant improvements in
larger, multi-block runs. The GPU also allows for intra-block parallelism,
dependant on the number of threads assigned to a block. Fig. 4.D shows how the
number of threads correlates with the simulation time of 82000 proposals of
Alt-MH for T0955. Finally, both forms of parallelism lead to lower overall
simulation time for the GPU implementation compared to the CPU implementation,
shown in Fig. 4.E. The CPU implementation (Intel I3-10100) was run for 1e5
proposals. The average time over 10 runs is plotted and the average minimum
energy was noted. The GPU was then run until a structure with the same minimum
energy is found for each protein.

This approach is intended to be similar to the time to solution metric
described in other literature [22, 15, 21]; here, there is an additional
challenge that the true ground state cannot be reached in reasonable time (or
likely ever due to error in the energy function). Hence we match GPU to CPU
optimality, rather than matching both to an optimality gap. Overall we observe
improvements in time to solution of 7.5x to 26.5. We predict these times can
be further improved with a custom hardware that caters more specifically to
MCMC processes than a GPGPU.

Figure 5:  Alt-MH is able to find low energy states at exceedingly large
$\beta$ for T0955 from CASP 13. $\beta$ is inverse temperature. Energy is
measured according to the distogram-only energy function. Each point is the
average of the minimum energies of 25 runs, where each run is $82000$
proposals. Shaded bands cover the range between the minimum and maximum values
found over 25 runs for each $\beta$.

A crucial benefit of Alternating Metropolis-Hastings is its’ diminished
dependence on temperature ($1/\beta$) when searching for low energy states. We
compare the performance of Alt-MH with traditional MH for varying $\beta$ in
Fig. 5. Localized MH prevents Alt-MH from locking into a sub-optimal state at
higher $\beta$. This is at odds with the behaviour of traditional Metropolis-
Hastings (MH), which can only reach (semi-)low energy states at a specific
$\beta$, ex. $\beta\approx 0.01$ in Fig. 5, above or below which the algorithm
fails. Alt-MH helps overcome the critical challenge of finding an ideal
$\beta$ or annealing schedule for new inputs.

##  IV Discussion

###  IV.1  Thoughts on Alternating Metropolis-Hastings

For Monte Carlo or Markov Chain Monte Carlo methods, visiting states with high
internal variation during a random walk generally requires proposal schemes
that can generate very specific, low probability states or running for
extended periods at a high temperature. An alternative approach is to allow
the algorithm to accept states that may increase global energy temporarily but
lead to accurate folding of local structures.

Localized MH fills an ideological role similar to that of noisy restarts in
gradient descent-based approaches. When the system converges to some local
minimum, it is necessary to introduce noise in such a way that the energy of
the system increases temporarily, so that it may once again minimize.
Localized MH acts as an intelligent restart, rather than a noisy restart, by
adding noise in a way that improves local structures. This makes it easier for
future iterations of the minimization scheme to find good global ground
states.

A critical benefit of the Alternating-MH algorithm is its diminished
dependence on temperature. Fig 5 shows that Alt-MH finds low energy states at
$\beta$s that far exceed the limits for traditional MH.

A point to note is that the existing standard for CPU protein simulations in
python, pyrosetta, is not built to efficiently find the energy of a subset of
a protein structure. Rather the full energy must be calculated, and the subset
must be gathered in post. Straightforward changes to the code base could allow
for significant improvements in the time required for alternating MH to
execute. If Alt-MH were made to run on par with the efficiency that gradient
descent based schemes do, then Time To Solution (TTS) measurements could be
made on CPU against gradient descent based schemes for the complex, refined
energy functions that pyrosetta allows for (trRosetta2).

###  IV.2  Other Monte Carlo Approaches to Protein Folding

Other Monte Carlo approaches to protein folding include population-based
schemes like EMC [23], genetic algorithms [24], fragment-swapping [25, 26,
14], rotamer substitutions [27], simulated annealing [28], and multi-particle
conformation sampling [29]. These methods have largely not been compared to or
applied on today’s state of the art protein folding pipelines. Recent works
have introduced efficient proposal schemes [30] that may integrate well with
the Alternating MH approach presented in this work.

Traditional Metropolis-Hastings has been used on all-atom force field
equations to sample ground states in a similar manner to this work [28]. This
work also begins at an unfolded protein and creates proposals by twisting
dihedral angles of residues in the chain. However, this is not applied to a
contact-map based energy function but rather as an alternative to molecular
dynamics on a more traditional inter-molecular force field energy function.
The idea of localized MH is also absent. Furthermore, this method is heavily
temperature dependent, while Alternating MH finds very low energy states
without requiring an ideal temperature, or an annealing schedule.

###  IV.3  Future Work

Incorporation of fragment libraries, novel proposal schemes, annealing
schedules, and further testing on parameters like the localization radius $d$
may yield improvements in minimization rate and effectiveness of the
Alternating MH algorithm. The Alternating MH algorithm may be applied to other
optimization problems like MaxCut or Travelling Salesman Problem (TSP); the
key will be to define what a localized segment means in the context of those
problems.

While we show that GPUs are able to accelerate Alt-MH over CPU, we believe
even greater improvements in time to solution could be possible with a
modified GPGPU that is catered specifically for Monte Carlo and Markov chain
Monte Carlo algorithms. Firstly, incorporating some number of local,
controllable RNG in each streaming multiprocessor could accelerate the drawing
of random numbers. Within the past few years, p-bits built out of stochastic
magnetic tunnel junctions have emerged as highly efficient, controllable
random number generators that could fill this need [31, 32]. We mention the
term ’controllable’ to motivate interconnected networks of p-bits which may
produce random numbers from arbitrary distributions. Secondly, limited shared
memory resources hinder block-wise MCMC parallelism. A trade-off of global and
local memory for increased shared memory resources would benefit the parallel
operation of MCMC algorithms where each block represents an MCMC chain. We
believe that these changes would also benefit the implementation of MCMC
algorithms for applications extending beyond protein structure prediction.

##  V  Conclusions

We present a novel Markov Chain Monte Carlo method termed Alternating
Metropolis-Hastings for optimizing energy functions, applied to protein
structure prediction. This algorithm is shown to outperform traditional MH,
and traditional MH with noisy restarts handily for 9 proteins taken from
recent CASP competitions. We also show that Alt-MH exhibits reduced dependence
on temperature compared with traditional MH, a crucial benefit when
considering novel problems for which an ideal temperature or schedule is
unknown. The algorithm is mapped to a GPGPU to show a 7.5x to 26.5x speedup.
Finally, we motivate the development of a modified GPGPU architecture with
per-block RNG and increased shared memory resources for accelerated MCMC
simulations.

##  Acknowledgements

The authors are grateful for many discussions with Shuvro Chowdhury, Jan
Kaiser, and Behtash Behin-Aein over the past year that played a key role in
shaping this work. The authors would also like to thank Kerem Camsari and Luke
Theogarajan for their advice in the early stages of this work. This work was
supported in part by ASCENT, one of six centers in JUMP, a Semiconductor
Research Corporation (SRC) program sponsored by DARPA.

## References

  * Anfinsen [1973] C. B. Anfinsen, Principles that govern the folding of protein chains, Science 181, 223 (1973).
  * Senior  _et al._ [2020] A. W. Senior, R. Evans, J. Jumper, J. Kirkpatrick, L. Sifre, T. Green, C. Qin, A. Žídek, A. W. R. Nelson, A. Bridgland, H. Penedones, S. Petersen, K. Simonyan, S. Crossan, P. Kohli, D. T. Jones, D. Silver, K. Kavukcuoglu, and D. Hassabis, Improved protein structure prediction using potentials from deep learning, Nature 577, 706 (2020), 1362 citations (Semantic Scholar/DOI) [2022-07-08].
  * Yang  _et al._ [2020] J. Yang, I. Anishchenko, H. Park, Z. Peng, S. Ovchinnikov, and D. Baker, Improved protein structure prediction using predicted interresidue orientations, Proceedings of the National Academy of Sciences 117, 1496 (2020), 635 citations (Semantic Scholar/DOI) [2022-07-08].
  * Marks  _et al._ [2011] D. S. Marks, L. J. Colwell, R. Sheridan, T. A. Hopf, A. Pagnani, R. Zecchina, and C. Sander, Protein 3D Structure Computed from Evolutionary Sequence Variation, PLoS ONE 6, e28766 (2011), 935 citations (Semantic Scholar/DOI) [2022-07-08].
  * Hopf  _et al._ [2014] T. A. Hopf, C. P. I. Schärfe, J. P. G. L. M. Rodrigues, A. G. Green, O. Kohlbacher, C. Sander, A. M. J. J. Bonvin, and D. S. Marks, Sequence co-evolution gives 3D contacts and structures of protein complexes, eLife 3, e03430 (2014), 394 citations (Semantic Scholar/DOI) [2022-07-08].
  * Di Lena  _et al._ [2012] P. Di Lena, K. Nagata, and P. Baldi, Deep architectures for protein contact map prediction, Bioinformatics 28, 2449 (2012), 263 citations (Semantic Scholar/DOI) [2022-07-08].
  * Li  _et al._ [2019] Y. Li, J. Hu, C. Zhang, D.-J. Yu, and Y. Zhang, ResPRE: high-accuracy protein contact prediction by coupling precision matrix with deep residual neural networks, Bioinformatics 35, 4647 (2019), 107 citations (Semantic Scholar/DOI) [2022-07-08].
  * Zheng  _et al._ [2019] W. Zheng, Q. Wuyun, Y. Li, S. M. Mortuza, C. Zhang, R. Pearce, J. Ruan, and Y. Zhang, Detecting distant-homology protein structures by aligning deep neural-network based contact maps, PLOS Computational Biology 15, e1007411 (2019), 31 citations (Semantic Scholar/DOI) [2022-07-08].
  * AlQuraishi [2019] M. AlQuraishi, AlphaFold at CASP13, Bioinformatics 35, 4862 (2019), 138 citations (Semantic Scholar/DOI) [2022-07-08].
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, On the limited memory BFGS method for large scale optimization, Mathematical Programming 45, 503 (1989), 5871 citations (Semantic Scholar/DOI) [2022-07-08].
  * Ma  _et al._ [2019] Y.-A. Ma, Y. Chen, C. Jin, N. Flammarion, and M. I. Jordan, Sampling can be faster than optimization, Proceedings of the National Academy of Sciences 116, 20881 (2019), 82 citations (Semantic Scholar/DOI) [2022-07-08].
  * Jumper  _et al._ [2021] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, A. Bridgland, C. Meyer, S. A. A. Kohl, A. J. Ballard, A. Cowie, B. Romera-Paredes, S. Nikolov, R. Jain, J. Adler, T. Back, S. Petersen, D. Reiman, E. Clancy, M. Zielinski, M. Steinegger, M. Pacholska, T. Berghammer, S. Bodenstein, D. Silver, O. Vinyals, A. W. Senior, K. Kavukcuoglu, P. Kohli, and D. Hassabis, Highly accurate protein structure prediction with AlphaFold, Nature 596, 583 (2021), 3193 citations (Semantic Scholar/DOI) [2022-07-08].
  * Baek  _et al._ [2021] M. Baek, F. DiMaio, I. Anishchenko, J. Dauparas, S. Ovchinnikov, G. R. Lee, J. Wang, Q. Cong, L. N. Kinch, R. D. Schaeffer, C. Millán, H. Park, C. Adams, C. R. Glassman, A. DeGiovanni, J. H. Pereira, A. V. Rodrigues, A. A. van Dijk, A. C. Ebrecht, D. J. Opperman, T. Sagmeister, C. Buhlheller, T. Pavkov-Keller, M. K. Rathinaswamy, U. Dalwadi, C. K. Yip, J. E. Burke, K. C. Garcia, N. V. Grishin, P. D. Adams, R. J. Read, and D. Baker, Accurate prediction of protein structures and interactions using a three-track neural network, Science 373, 871 (2021), 505 citations (Semantic Scholar/DOI) [2022-07-08].
  * Senior  _et al._ [2019] A. W. Senior, R. Evans, J. Jumper, J. Kirkpatrick, L. Sifre, T. Green, C. Qin, A. Žídek, A. W. R. Nelson, A. Bridgland, H. Penedones, S. Petersen, K. Simonyan, S. Crossan, P. Kohli, D. T. Jones, D. Silver, K. Kavukcuoglu, and D. Hassabis, Protein structure prediction using multiple deep neural networks in the 13th Critical Assessment of Protein Structure Prediction (CASP13), Proteins: Structure, Function, and Bioinformatics 87, 1141 (2019), 148 citations (Semantic Scholar/DOI) [2022-07-08].
  * Kaiser  _et al._ [2021] J. Kaiser, R. Jaiswal, B. Behin-Aein, and S. Datta, Benchmarking a Probabilistic Coprocessor (2021), 3 citations (Semantic Scholar/arXiv) [2022-07-08] arXiv:2109.14801 [cond-mat].
  * Hansmann [1997] U. H. E. Hansmann, Parallel tempering algorithm for conformational studies of biological molecules, Chemical Physics Letters 281, 140 (1997), 903 citations (Semantic Scholar/DOI) [2022-07-29].
  * Liu  _et al._ [2000] J. S. Liu, F. Liang, and W. H. Wong, The Multiple-Try Method and Local Optimization in Metropolis Sampling, Journal of the American Statistical Association 95, 121 (2000), publisher: [American Statistical Association, Taylor & Francis, Ltd.].
  * Liang and Wong [2000] F. Liang and W. H. Wong, EVOLUTIONARY MONTE CARLO: APPLICATIONS TO C p MODEL SAMPLING AND CHANGE POINT PROBLEM, Statistica Sinica 10, 317 (2000), publisher: Institute of Statistical Science, Academia Sinica.
  * Craiu  _et al._ [2009] R. V. Craiu, J. Rosenthal, and C. Yang, Learn From Thy Neighbor: Parallel-Chain and Regional Adaptive MCMC, Journal of the American Statistical Association 104, 1454 (2009), publisher: [American Statistical Association, Taylor & Francis, Ltd.].
  * Grenander and Miller [1994] U. Grenander and M. I. Miller, Representations of Knowledge in Complex Systems, Journal of the Royal Statistical Society. Series B (Methodological) 56, 549 (1994), publisher: [Royal Statistical Society, Wiley].
  * Kaiser and Datta [2021] J. Kaiser and S. Datta, Probabilistic computing with p-bits, Applied Physics Letters 119, 150503 (2021), 0 citations (Semantic Scholar/DOI) [2022-07-08].
  * Patel  _et al._ [2020] S. Patel, L. Chen, P. Canoza, and S. Salahuddin, Ising Model Optimization Problems on a FPGA Accelerated Restricted Boltzmann Machine (2020), 14 citations (Semantic Scholar/arXiv) [2022-07-30] arXiv:2008.04436 [physics].
  * Liang and Wong [2001] F. Liang and W. H. Wong, Evolutionary Monte Carlo for protein folding simulations, The Journal of Chemical Physics 115, 3374 (2001), 187 citations (Semantic Scholar/DOI) [2022-07-08].
  * Hansmann and Okamoto [1999] U. H. Hansmann and Y. Okamoto, New Monte Carlo algorithms for protein folding, Current Opinion in Structural Biology 9, 177 (1999), 271 citations (Semantic Scholar/DOI) [2022-07-08].
  * Park  _et al._ [2018] H. Park, S. Ovchinnikov, D. E. Kim, F. DiMaio, and D. Baker, Protein homology model refinement by large-scale energy optimization, Proceedings of the National Academy of Sciences 115, 3054 (2018), 52 citations (Semantic Scholar/DOI) [2022-07-08].
  * Kandathil  _et al._ [2018] S. M. Kandathil, M. Garza-Fabre, J. Handl, and S. C. Lovell, Improved fragment-based protein structure prediction by redesign of search heuristics, Scientific Reports 8, 13694 (2018), 9 citations (Semantic Scholar/DOI) [2022-07-08].
  * Kuhlman and Bradley [2019] B. Kuhlman and P. Bradley, Advances in protein structure prediction and design, Nature Reviews Molecular Cell Biology 20, 681 (2019), 263 citations (Semantic Scholar/DOI) [2022-07-08].
  * Heilmann  _et al._ [2020] N. Heilmann, M. Wolf, M. Kozlowska, E. Sedghamiz, J. Setzler, M. Brieg, and W. Wenzel, Sampling of the conformational landscape of small proteins with Monte Carlo methods, Scientific Reports 10, 18211 (2020), 5 citations (Semantic Scholar/DOI) [2022-07-08].
  * Wong  _et al._ [2018] S. W. K. Wong, J. S. Liu, and S. C. Kou, Exploring the conformational space for protein folding with sequential Monte Carlo, The Annals of Applied Statistics 12, 10.1214/17-AOAS1124 (2018), 5 citations (Semantic Scholar/DOI) [2022-07-08].
  * Cabeza de Vaca  _et al._ [2018] I. Cabeza de Vaca, Y. Qian, J. Z. Vilseck, J. Tirado-Rives, and W. L. Jorgensen, Enhanced Monte Carlo Methods for Modeling Proteins Including Computation of Absolute Free Energies of Binding, Journal of Chemical Theory and Computation 14, 3279 (2018), 21 citations (Semantic Scholar/DOI) [2022-07-08].
  * Camsari  _et al._ [2017] K. Y. Camsari, R. Faria, B. M. Sutton, and S. Datta, Stochastic p -Bits for Invertible Logic, Physical Review X 7, 031014 (2017), 156 citations (Semantic Scholar/DOI) [2022-08-02].
  * Camsari  _et al._ [2019] K. Y. Camsari, B. M. Sutton, and S. Datta, p-bits for probabilistic spin logic, Applied Physics Reviews 6, 011305 (2019), 55 citations (Semantic Scholar/DOI) [2022-08-02].