11institutetext:  Johns Hopkins University, Baltimore, MD, 21218, USA
22institutetext: University of Maryland School of Medicine and R Adams Cowley
Shock Trauma Center, Baltimore, MD, 21201, USA
22email<EMAIL_ADDRESS>
# FastSAM-3DSlicer: A 3D-Slicer Extension for 3D Volumetric Segment Anything
Model with Uncertainty Quantification

Yiqing Shen 11 Xinyuan Shao 11 Blanca Inigo Romillo 11 David Dreizin 22
Mathias Unberath(✉) 11

###### Abstract

Accurate segmentation of anatomical structures and pathological regions in
medical images is crucial for diagnosis, treatment planning, and disease
monitoring. While the Segment Anything Model (SAM) and its variants have
demonstrated impressive interactive segmentation capabilities on image types
not seen during training without the need for domain adaptation or retraining,
their practical application in volumetric 3D medical imaging workflows has
been hindered by the lack of a user-friendly interface. To address this
challenge, we introduce FastSAM-3DSlicer, a 3D Slicer extension that
integrates both 2D and 3D SAM models, including SAM-Med2D, MedSAM, SAM-Med3D,
and FastSAM-3D. Building on the well-established open-source 3D Slicer
platform, our extension enables efficient, real-time segmentation of 3D
volumetric medical images, with seamless interaction and visualization. By
automating the handling of raw image data, user prompts, and segmented masks,
FastSAM-3DSlicer provides a streamlined, user-friendly interface that can be
easily incorporated into medical image analysis workflows. Performance
evaluations reveal that the FastSAM-3DSlicer extension running FastSAM-3D
achieves low inference times of only 1.09 seconds per volume on CPU and 0.73
seconds per volume on GPU, making it well-suited for real-time interactive
segmentation. Moreover, we introduce an uncertainty quantification scheme that
leverages the rapid inference capabilities of FastSAM-3D for practical
implementation, further enhancing its reliability and applicability in medical
settings. FastSAM-3DSlicer offers an interactive platform and user interface
for 2D and 3D interactive volumetric medical image segmentation, offering a
powerful combination of efficiency, precision, and ease of use with SAMs. The
source code and a video demonstration are publicly available at
https://github.com/arcadelab/FastSAM3D_slicer.

###### Keywords:

Foundation Model Deep Learning Segment Anything Model (SAM) Interactive
Segmentation Interface 3D Slicer.

##  1 Introduction

Precise segmentation of anatomical structures and pathological regions from
medical images is essential for accurate diagnosis, treatment planning, and
monitoring of disease progression [4]. However, manual segmentation is time-
consuming, labor-intensive, and prone to inter-observer variability. Deep-
learning-driven automatic segmentation models have shown promise in reducing
manual effort, they often struggle to generalize across diverse datasets,
anatomical variations, and unseen pathologies during inference [6]. The
Segment Anything Model (SAM) [5] and its variants for volumetric medical
images, such as SAM-Med3D [11], have emerged as flexible zero-shot solutions
that enable interactive segmentation of novel objects without requiring
retraining. These models are designed to generalize across various tasks and
datasets through large-scale pre-training and support for manual prompting.
The interactive nature of SAM-based segmentation makes fast inference times
particularly important, as it allows users to make immediate adjustments and
corrections, improving the accuracy and efficiency of the segmentation
process. FastSAM-3D [10], a computationally efficient 3D SAM architecture, has
been specifically optimized for real-time interactive segmentation of 3D
volumes such as computed tomography (CT). FastSAM-3D utilizes a compact Vision
Transformer (ViT) encoder [2], distilled from the larger SAM-Med3D [11], and
incorporates an efficient 3D Sparse Flash Attention [10] mechanism to reduce
computational costs while maintaining high segmentation quality.

Despite the potential of efficient SAM variants to enable highly responsive
interactive segmentation, there is currently a lack of user-friendly
interfaces to facilitate their practical application in medical image analysis
workflow. Interactive prompting in 3D medical volumes is more challenging
compared to 2D interfaces for natural images, due to the increased complexity
of visualizing 3D data on a 2D screen and the need for seamless integration
with existing medical image analysis workflows. 3D
Slicer111https://www.slicer.org [3] is an open-source software platform widely
used for the analysis and visualization of volumetric medical images. It
supports a variety of plug-ins and offers a robust framework for the
integration of new tools, making it an ideal choice for creating an interface
for interactive volumetric image segmentation, e. g., with FastSAM-3D. In this
manuscript, we describe FastSAM-3DSlicer, a plugin for integrating 2D and 3D
SAM models, including the efficient FastSAM-3D, into the well-established
image analysis platform 3D Slicer. This extension enables users to load 3D
volumes, interactively annotate structures of interest using point prompts,
and visualize the resulting segmentations in real time within a familiar
software environment.

In summary, our contributions are three-fold. Firstly, we propose a novel 3D
Slicer-based extension for both 2D and 3D SAM models, including the efficient
FastSAM-3D, for volumetric image segmentation. It enables interactive
prompting in a 3D manner with SAM. Secondly, we show the quantitative
comparison of the inference time of different SAM models within our interface
on both CPU and GPU environments. Finally, we propose an uncertainty
quantification scheme based on the fast inference speed of FastSAM-3D in our
extension. It can guide the user for better prompting.

Figure 1:  An illustration of the interface of our 3D Slicer extension i.e.,
FastSAM-3DSlicer. It allows users to load 3D volumes (red), select the mode
for uncertainty quantification (green), choose SAM models (brown), manage
masks (blue) and prompts (yellow), and visualize segmentation and input volume
in real-time (purple). This integration enhances image analysis workflows by
providing an intuitive and efficient interface for the interactive annotation
and segmentation of volumetric medical images with SAM.  Figure 2:  The
overall implementation workflow of FastSAM-3DSlicer. When FastSAM-3DSlicer is
started for the first time, it runs the configure dependency block to check
and install all required dependencies. The setup block runs every time
FastSAM-3DSlicer is started to initialize parameters, models, and the user
interface. Users can select a SAM model using a radio button. The green blocks
represent the components of the 3D Slicer environment, including the volume
node, markup nodes, and segmentation node, which handle the raw image data,
user interactions, and generated masks, respectively. The blue blocks
represent the FastSAM-3DSlicer Python components, which include configuring
dependencies, setup, and the various supported 2D and 3D SAM models (SAM-Med2D
[1], MedSAM [8], SAM-Med3D [11], and FastSAM-3D [10]). Black text indicates
processes or actions within the workflow, while gray text represents data
types or storage formats used throughout the process.

##  2 Methods

###  2.0.1 Overview of the 3D Slicer Extension

The interface for our FastSAM-3DSlicer extension is illustrated in Fig. 1,
with its overall implementation details shown in Fig. 2. Unlike previous works
integrating 2D SAMs into 3D Slicer, such as TomoSAM [9] and SAMME [7],
FastSAM-3D’s efficiency at the volume level eliminates the need to prepare
image embeddings before interactive segmentation (i.e., user prompting) in 3D
Slicer. This improvement enhances the user experience by reducing pre-
processing time and allowing for more immediate interaction with the data.
Upon importing a 3D volumetric image into FastSAM-3DSlicer, the extension
automatically generates three node types, including a volume node, a
segmentation node, and two markup nodes. These node types represent specific
data structures within 3D Slicer that manage different aspects of the image
and segmentation process. The volume node contains the raw 3D image in
grayscale, represented as a NumPy array, which serves as the input for the SAM
model The segmentation node stores all the segmented masks, which are updated
when users add new point prompts and generated when users perform the add mask
operation. The two markup nodes hold all the include and exclude points that
the user can add interactively. When a user adds an include or exclude point
prompt, FastSAM-3DSlicer first converts this point from RAS to XYZ coordinates
based on the affine matrix stored in the input NIfTI file. Using this new
point prompt and previous input point prompts, FastSAM-3DSlicer crops or pads
the raw image to match the selected SAM model’s input size and feeds it into
the image encoder to generate image embeddings in real time. The prompt
encoder processes all the input points to generate prompt embeddings. The mask
decoder then translates the image embeddings, prompt embeddings, and previous
masks into a new mask. This segmentation mask is resized to the raw image
dimensions based on saved coordinate information and updated in the
segmentation node. In addition to supporting the FastSAM-3D model,
FastSAM-3DSlicer also supports SAM-Med3D [11], MedSAM [8], and SAM-Med2D [1],
all of which follow the same structural process.

###  2.0.2 Uncertainty Quantification Scheme

Uncertainty quantification in FastSAM-3DSlicer provides users with a measure
of confidence in the segmentation results, which can be used to guide user
prompting. Regions with higher uncertainty indicate a greater need for
additional prompts. Our method leverages the efficiency of FastSAM-3D by
running the image encoder once and performing multiple decoding steps to
generate an ensemble of segmentations. It begins with the initial
segmentation, where the image encoder processes the input 3D volumetric image
$\mathbf{I}$ to produce its image embedding
$\mathbf{E}=\text{Encoder}(\mathbf{I})$. Next, the mask decoder generates the
initial segmentation logits $\widehat{\mathbf{M}}_{0}$ based on the image
embedding $\mathbf{E}$ and the initial set of point prompts ($\mathbf{P}_{0}$)
provided from the user, i.e.
$\widehat{\mathbf{M}}_{0}=\text{Decoder}(\mathbf{E},\mathbf{P}_{0})$. The
segmentation mask $\mathbf{M}_{0}$ is the binarized segmentation logits
$\widehat{\mathbf{M}}_{0}$, obtained by applying a threshold $\tau$:

$\mathbf{M}_{0}=\mathbbm{1}(\widehat{\mathbf{M}}_{0}>\tau),$ (1)

where $\mathbbm{1}(\cdot)$ is the indicator function. To quantify uncertainty,
subsequent point prompts are sampled from the initial segmentation mask
$\mathbf{M}_{0}$. These pseudo-point prompts are used to run the decoder
multiple times, each time producing a slightly different segmentation mask due
to variations in the sampled prompts. Let $\mathbf{P}_{i}$ denote the point
prompts sampled from the initial segmentation mask $\mathbf{M}_{0}$ for the
$i$th iteration: $\mathbf{P}_{i}=\text{SamplePrompts}(\mathbf{M}_{0})$. The
decoder is then run $N$ times using these sampled prompts, while keeping the
image encoder constant, to produce $N$ different segmentation masks
$\\{\widehat{\mathbf{M}}_{i}\\}_{i=1}^{N}$ with
$\widehat{\mathbf{M}}_{i}=\text{Decoder}(\mathbf{E},\mathbf{P}_{i})$. Since
the image encoder only runs once, the majority of the computational efficiency
is preserved, as the decoder accounts for a smaller proportion of the total
computation. Inspired by the self-ensembling [13], the segmentation logits
from each decoder run are averaged to produce the final ensemble result
$\widehat{\mathbf{M}}$:

$\widehat{\mathbf{M}}=\frac{1}{N}\sum_{i=1}^{N}\widehat{\mathbf{M}}_{i}.$ (2)

This averaging process not only provides a robust final segmentation mask but
also allows for the calculation of uncertainty. Specifically, the variability
among the $N$ segmentation masks can be used to estimate the uncertainty,
formally expressed as the standard deviation or variance of the logits at each
voxel:

$\text{Uncertainty}(x)=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(\widehat{\mathbf{M}}_{i}(x)-\widehat{\mathbf{M}}(x))^{2}}$
(3)

where $x$ denotes a voxel in the 3D volume.

##  3 Experiments

###  3.0.1 Implementation Details

We implemented the proposed 3D Slicer extension using 3D Slicer version 5.6.2
and Python version 3.10. The experiments were conducted in two distinct
environments to evaluate inference times. The first environment was a CPU-only
setup utilizing a laptop-level AMD Ryzen 5 5500U CPU, while the second
environment employed a GPU setup with one NVIDIA RTX 2060 GPU. Following
previous work [10], we use the test set of totalsegmentator [12] to test the
inference time. All data are prepared in NIfTI format for loading. Code for
the 3D Slicer extension is available at
https://github.com/arcadelab/FastSAM3D_slicer. A video demo for our 3D slicer
extension is available at https://www.youtube.com/watch?v=oJ9ZhnPWqSs.

###  3.0.2 Results for Inference Time

Table 1: Comparison of inference times for different SAMs on CPU and GPU in our 3D Slicer extension. Results in red indicate slice-level inference times, while results in blue indicate volume-level inference times.  Models |  Inference Time w/ CPU |  Inference Time w/ GPU
---|---|---
SAM-Med2D [1] |  1.52 seconds per slice |  0.52 seconds per slice
MedSAM [8] |  48.9 seconds per slice |  12.69 seconds per slice
SAM-Med3D [11] |  7.75 seconds per volume |  1.76 seconds per volume
FastSAM-3D [10] |  1.09 seconds per volume |  0.73 seconds per volume

Table 1 presents a comparison of inference times for different SAMs with both
CPU and GPU environments within our FastSAM-3DSlicer extension. The SAM-Med2D
[1] exhibits an inference time of 1.52 seconds per slice on the CPU and 0.52
seconds per slice on the GPU. While SAM-Med2D [1] is effective for 2D slice-
based segmentation, it fails to address volume-level segmentation in 3D
Slicer, as the user needs to provide prompts for each slice individually,
which results in longer inference time. MedSAM [8] shows considerably higher
inference times, with 48.9 seconds per slice on the CPU and 12.69 seconds per
slice on the GPU. The increased processing time is due to its higher
resolution of $1024\times 1024$ compared to the $256\times 256$ resolution of
SAM-Med2D [1]. Its high computational cost limits its practicality for real-
time applications in image analysis settings. The SAM-Med3D [11], optimized
for 3D segmentation, achieves 7.75 seconds per volume on the CPU and 1.76
seconds per volume on the GPU. This demonstrates a substantial improvement
over MedSAM [8], particularly in GPU environments, making it a more viable
option for volume-level segmentation. FastSAM-3D [10] significantly
outperforms the other models in terms of inference speed ($p<0.01$). It
achieves 1.09 seconds per volume on the CPU and 0.73 seconds per volume on the
GPU. This reduction in inference time highlights the efficiency and
optimization of FastSAM-3D for real-time interactive segmentation of 3D
volumes. FastSAM-3D demonstrates the lowest inference times on volume-level
segmentation across both CPU and GPU environments, affirming its suitability
for integration into image analysis workflows where speed and efficiency are
critical. The improvements in inference time not only facilitate faster
segmentation but also enable more immediate and iterative interaction with the
volumetric data, thereby enhancing the overall utility of the 3D Slicer
extension in medical image analysis applications.

###  3.0.3 Illustrative Visual Examples

Figure 3: Qualitative comparison of segmentation results from different SAM
models in the FastSAM-3DSlicer interface. Each row showcases the segmentation
performance on a specific anatomical structure or region of interest. The red
masks represent the ground truth segmentations, while the green masks depict
the predicted segmentations from each model.

Fig. 3 presents a qualitative comparison of segmentation results obtained from
different SAM models integrated into the FastSAM-3DSlicer interface. The
examples showcase the performance of each model on various anatomical
structures and regions of interest. FastSAM-3D and SAM-Med3D generate
segmentations for the entire 3D volume, demonstrating their ability to capture
spatial context and produce coherent masks. In contrast, MedSAM and SAM-Med2D
operate on 2D slices, resulting in smaller and more localized segmentations
when visualized in the 3D view. Across all examples, FastSAM-3D exhibits the
highest agreement with the ground truth, highlighting its superior performance
in terms of both efficiency and accuracy for real-time 3D interactive
segmentation within the FastSAM-3DSlicer.

###  3.0.4 Results for Uncertainty Quantification

Figure 4: Visualization of the uncertainty quantification results obtained
through self-ensembling in the FastSAM-3DSlicer. Each column represents a
different SAM model, while the last column shows the ground truth
segmentation. The uncertainty maps are displayed as heatmaps overlaid on the
original images, where lighter colors indicate higher uncertainty (lower
probability) in the segmentation.

Fig. 4 showcases the uncertainty quantification results obtained through the
self-ensembling approach in the FastSAM-3DSlicer extension. The results
demonstrate that FastSAM-3D exhibits the lowest overall uncertainty among the
compared models, highlighting its robustness and reliability for real-time 3D
interactive segmentation. It allows for the estimation of uncertainty by
calculating the variance among the ensemble of segmentations. The uncertainty
information provided by FastSAM-3DSlicer is particularly valuable for guiding
user interactions, as it identifies regions where additional user prompts may
be required to improve segmentation accuracy. By focusing on areas of high
uncertainty, users can iteratively refine the segmentation results, leading to
more precise and reliable delineations of anatomical structures or regions of
interest.

##  4 Conclusion

We presented FastSAM-3DSlicer, a 3D Slicer extension designed to facilitate
real-time, interactive segmentation of 3D volumetric medical images with SAMs.
Our extension integrates SAM-Med2D, MedSAM, SAM-Med3D, and FastSAM-3D,
offering a user-friendly interface that automates the handling of user prompts
and segmented masks within the familiar 3D Slicer environment. Our extension
demonstrates superior performance in terms of inference time, particularly
with the FastSAM-3D model, which achieves low inference times on both CPU and
GPU environments. This makes FastSAM-3D highly suitable for real-time
applications, reducing the computational burden while maintaining high
segmentation quality. Furthermore, the integration of an innovative
uncertainty quantification scheme leverages the rapid inference capabilities
of FastSAM-3D, providing users with additional information about the
reliability of the segmentation results. Overall, by combining computational
efficiency, precision, and ease of use, FastSAM-3DSlicer addresses the need
for a user-friendly interface for SAMs, thereby enhancing decision-making
processes and improving patient outcomes in medical settings.

## References

  * [1] Cheng, J., et al.: Sam-med2d. arXiv preprint arXiv:2308.16184 (2023)
  * [2] Dosovitskiy, A., et al.: An image is worth 16x16 words: Transformers for image recognition at scale. arXiv preprint arXiv:2010.11929 (2020)
  * [3] Fedorov, A., et al.: 3d slicer as an image computing platform for the quantitative imaging network. Magnetic resonance imaging 30(9), 1323–1341 (2012)
  * [4] Isensee, F., et al.: nnu-net: a self-configuring method for deep learning-based biomedical image segmentation. Nature methods 18(2), 203–211 (2021)
  * [5] Kirillov, A., et al.: Segment anything. arXiv preprint arXiv:2304.02643 (2023)
  * [6] Liu, X., et al.: Review of deep learning based automatic segmentation for lung cancer radiotherapy. Frontiers in oncology 11, 717039 (2021)
  * [7] Liu, Y., et al.: Segment any medical model extended. In: Medical Imaging 2024: Image Processing. vol. 12926, pp. 411–422. SPIE (2024)
  * [8] Ma, J., et al.: Segment anything in medical images. Nature Communications 15(1), 654 (2024)
  * [9] Semeraro, F., et al.: Tomosam: a 3d slicer extension using sam for tomography segmentation. arXiv preprint arXiv:2306.08609 (2023)
  * [10] Shen, Y., et al.: Fastsam3d: An efficient segment anything model for 3d volumetric medical images. arXiv preprint arXiv:2403.09827 (2024)
  * [11] Wang, H., et al.: Sam-med3d. arXiv preprint arXiv:2310.15161 (2023)
  * [12] Wasserthal, J., et al.: Totalsegmentator: Robust segmentation of 104 anatomic structures in ct images. Radiology: Artificial Intelligence 5(5) (2023)
  * [13] Zhang, Y., Hu, S., Jiang, C., Cheng, Y., Qi, Y.: Segment anything model with uncertainty rectification for auto-prompting medical image segmentation. arXiv preprint arXiv:2311.10529 (2023)


10

#  PRESENT STATUS OF PICOLON PROJECT

K. Fushimi∗ Division of Science and Technology, Tokushima University, 2-1
Minami Josanajima-cho, Tokushima city, Tokushima , 770-8506, Japan  D.
Chernyak  Department of Physics and Astronomy, University of Alabama,
Tuscaloosa, Alabama 35487, USA  H. Ejiri  Research Center for Nuclear Physics,
Osaka University, 10-1 Mihogaoka Ibaraki city, Osaka, 567-0042, Japan  K. Hata
Research Center for Neutrino Science, Tohoku University, 6-3 Aramaki Aza Aoba,
Aobaku, Sendai city, Miyagi, 980-8578, Japan  R. Hazama  Department of
Environmental Science and Technology, Osaka Sangyo University, 3-1-1
Nakagaito, Daito city, Osaka, 574-8530, Japan  T. Iida  Faculty of Pure and
Applied Sciences, University of Tsukuba, 1-1-1 Tennoudai, Tsukuba city,
Ibaraki, 305-8571, Japan  H. Ikeda  K. Imagawa  I. S. C. Lab. , 5-15-24
Torikai Honmachi, Settsu city, Osaka, 566-0052, Japan  K. Inoue  H. Ito
Department of Physics, Faculty of Science and Technology, Tokyo University of
Science, Noda, Chiba 278-8510, Japan  T. Kisimoto  Department of Physics,
Osaka University, 1-1 Machikaneyama-cho, Toyonaka city, Osaka 560-0043, Japan
M. Koga  K. Kotera  Graduate School of Advanced Technology and Science,
Tokushima University, 2-1 Minami Josanajima-cho, Tokushima city, Tokushima ,
770-8506, Japan  A. Kozlov  National Research Nuclear University “MEPhI”
(Moscow Engineering Physics Institute), Moscow, 115409, Russia  S. Kurosawa
Institute for Materials Research, Tohoku University, 2-1-1 Katahira, Aoba-ku,
Sendai, Miyagi 980-8577, Japan  Institute of Laser Engineering, Osaka
University, 2-6 Yamadaoka, Suita, Osaka 565-0871, Japan  K. Nakamura  Butsuryo
College of Osaka, 3-33 Ohtori Kitamachi, Nishi ward, Sakai city, Osaka,
593-8328, Japan  R. Orito  A. Sakaguchi  A. Sakaue  T. Shima  Y. Takaku  Y.
Takemoto  Institute of Laser Engineering, Osaka University, 2-6 Yamadaoka,
Suita, Osaka 565-0871, Japan  S. Umehara  Y. Urano  Y. Yamamoto  K. Yasuda  S.
Yoshida

###### Abstract

The existence of cosmic dark matter and neutrino properties are long-standing
problems in cosmology and particle physics. These problems have been
investigated by using radiation detectors. We will discuss the application of
inorganic crystal scintillators to studies on dark matter and neutrino
properties. A large volume and high-purity inorganic crystal is a promising
detector for investigating dark matter and neutrino.

##  1 Introduction

###  1.1 Dark matter problem

Cosmic dark matter is an unknown elementary particle that does not interact
via electromagnetic interaction. There are two attractive candidates from
particle physics. One is weakly interacting massive particles (WIMPs), and the
other is the axion.

WIMPs are a sort of hypothetical elementary particles which interact via weak
interaction. Their mass ranges from sub-GeV to a few TeV. The root mean square
of the WIMPs velocity is 238 km/s [1, 2, 3]. It gives recoil energy to an
atomic nucleus by elastic scattering. The mean recoil energy and event rate
are modulated annually because of the proper motion of the sun and the earth’s
revolution. The maximum of WIMPs energy is given at the beginning of June, and
the minimum is given at the beginning of December [4, 5]. Figure 1 shows the
expected energy spectrum of Na nuclei in a NaI(Tl) scintillator recoiled by
WIMPs with its mass of 50 GeV$/c^{2}$, where $c$ is the speed of light. The
energy spectra were calculated with the cross-section of elastic scattering
between proton and WIMPs as $\sigma_{\mathrm{p}-\chi}=10^{-36}$ cm2.

Figure 1: The expected energy spectra of recoil Na in NaI(Tl) scintillator.
The black and red lines are the recoil energy spectra observed in June and
December.

Axion is another candidate for dark matter. It was proposed to solve the CP
violation problem in strong interaction [6, 7]. The axion is a light boson and
interacts via the Peccei-Quinn mechanism, which is analogous to
photoelectronic interaction.

Many experimental works have been done to detect dark matter candidates by
using various radiation detectors. The semiconductor detector has a
significant advantage in particle identification. The ratio of phonons to
ionization is small for electrons and large for recoil nuclei[8]. The
semiconductor detectors can perform almost background-free experiments since
the WIMPs signal gives nuclear recoil signal.

The other detector, which gives excellent results, is a large-volume xenon
scintillator [9, 10]. The two-phase xenon detector not only discriminates
electron-recoil and nuclear recoil but also identifies the position of the
interaction. XENONnT detector combines a scintillation detector and a TPC
(time projection chamber) detector. TPC measures the ionization signal and the
position of interaction. The signal intensity ratio between ionization and
scintillation enables discriminating between nuclear and electron recoil
events. Recently, XENONnT group reported the lowest background level in
electron recoil events as $(4.3\pm 0.4)\times 10^{-5}$ day-1kg-1keV-1 by $5.9$
ton liquid xenon detector [10]. Despite the low-background detectors, no
significant signal of WIMPs has been reported. The experimental limit on the
cross-section was dug down to $10^{-45}$ cm2. We need to develop a low-
background detector as low as $10^{-5}$ day-1kg-1keV-1 .

###  1.2 Neutrino mass

The neutrino mass is not measured directly, although the evidence of a finite
mass was reported by neutrino oscillation. The neutrinos have three flavor
eigenstates in weak interaction, $\nu_{\mathrm{e}}$, $\nu_{\mu}$, and
$\nu_{\tau}$. If the neutrinos have finite mass, the flavor eigenstate is
expressed by a linear combination of the mass eigenstate, $\nu_{1}$,
$\nu_{2}$, and $\nu_{3}$. The relation between flavor eigenstate and mass
eigenstate is expressed as

$\begin{pmatrix}\nu_{\mathrm{e}}\\\ \nu_{\mu}\\\ \nu_{\tau}\\\
\end{pmatrix}=\begin{pmatrix}U_{\mathrm{e}1}&U_{\mathrm{e}2}&U_{\mathrm{e}3}\\\
U_{\mu 1}&U_{\mu 2}&U_{\mu 3}\\\ U_{\tau 1}&U_{\tau 2}&U_{\tau 3}\\\
\end{pmatrix}\begin{pmatrix}\nu_{1}\\\ \nu_{2}\\\ \nu_{3}\\\ \end{pmatrix}$
(1)

The observable eigenstate of the neutrino is its weak eigenstate because we
observe them via weak interaction with electrons and muons. The neutrino
flavor changes traveling between the emission point and detection point. The
probability of detecting an electron-neutrino at a distance depends on the
difference in squared mass

$\Delta m^{2}_{ij}=m^{2}_{i}-m^{2}_{j},$ (2)

where $m_{i},i=1,2,3$ is the mass of the neutrinos $\nu_{i},i=1,2,3$. The
detection probability of the original neutrino flavor oscillates with the
distance of the neutrino source, called neutrino oscillation. The neutrino
oscillation insists that the neutrino has a finite mass. However, the
experimental results do not tell its absolute value.

The neutrino-less double beta decay ($0\nu\beta\beta$) is one of the most
sensitive methods to measure the absolute mass of the neutrino [11, 12].
Double beta decay occurs in the nuclei if single beta decay is prohibited and
double beta decay is allowed. The nucleons that make up the nucleus are
arranged in order from the lowest energy orbital. Two nucleons in the same
orbit settle in a lower energy state when their spins are opposite. Nuclei
with even numbers of protons and neutrons, called even-even nuclei, have less
mass than odd-odd nuclei with an odd number of both. As a result, 136Cs next
to 136Xe is heavier than 136Xe, and 136Ba next to it is lighter than 136Xe.
This energy difference causes double beta decay.

The two-neutrino double beta decay ($2\nu\beta\beta$) is expressed by the
decay formula as

$(A,Z)\rightarrow(A,Z+2)+2\mathrm{e}^{-}+\bar{\nu_{\mathrm{e}}},$ (3)

where $A$ and $Z$ are the nucleus’s mass number and atomic number (Left of
Figure 2). The $2\nu\beta\beta$ occurs in the flame of the standard model of
particle physics.

An electron has two chirality states, right-handed ($\mathrm{e}_{\mathrm{R}}$)
and left-handed ($\mathrm{e}_{\mathrm{L}}$). Chirality is defined as

$\displaystyle\psi_{\mathrm{L}}$ $\displaystyle=\frac{1-\gamma^{5}}{2}\psi$
(4) $\displaystyle\psi_{\mathrm{R}}$
$\displaystyle=\frac{1+\gamma^{5}}{2}\psi.$ (5)

The four-component Dirac spinor $\psi$ is written in terms of two upper
components of the left chirality field $\psi_{\mathrm{L}}$ and two lower
components of the right chirality field $\psi_{\mathrm{R}}$. The charged
leptons have two chirality states, L and R. On the other hand, neutrinos have
only left-handed chirality. Therefore, only left-handed neutrino
$\nu_{\mathrm{L}}$ and right-handed antineutrino $\bar{\nu_{\mathrm{R}}}$
exist in the standard model.

Figure 2: Image of the double beta decay. The neutron and proton in a nucleus
are blue and red circles. A solid line and a dashed line stand for the
electron and neutrino, respectively

$0\nu\beta\beta$ process gives a significant information to the beyond
standard model. $0\nu\beta\beta$ process occurs if the neutrino has right-
handed chirality and the neutrino and anti-neutrino are the same, called
Majorana neutrino. The $0\nu\beta\beta$ is expressed as

$(A,Z)\rightarrow(A,Z+2)+2\mathrm{e}^{-}.$ (6)

In this case, a neutrino emitted by one neutron is absorbed by another in the
same nucleus or annihilates each other (see right side of Figure 2). The decay
probability of $0\nu\beta\beta$ is expressed by Majorana mass
$\left<m\right>$, and right-handed interaction in terms of
$\left<\lambda\right>$ and $\left<\eta\right>$. The half-life of
$0\nu\beta\beta$ decay, $T^{0\nu}_{1/2}$, is expressed as

$\frac{\ln 2}{T^{0\nu}_{1/2}}=G^{0\nu}\left|M^{0\nu}\right|^{2}K,$ (7)

where $G^{0\nu}$ is the kinematical factor and $M^{0\nu}$ is the nuclear
matrix element. The kinematical factor and nuclear matrix element are obtained
by nuclear physics [11, 12]. $K$ is the neutrino mass term
$\left<m_{\nu}\right>$, and the right-handed current terms
$\left<\lambda\right>$ and $\left<\eta\right>$.

The expected half-life of $0\nu\beta\beta$ is more than $10^{6}$ longer than
that of $2\nu\beta\beta$ process. KamLAND-Zen group searched for
$0\nu\beta\beta$ of 136Xe by 970 kg$\cdot$yr exposure [13]. They obtained a
lower limit on the half-life of $0\nu\beta\beta$ decay of
$T_{1/2}^{0\nu}>2.3\times 10^{26}$ yr at 90 % C.L. Their limit corresponds to
upper limits on the effective Majorana neutrino mass of $36\sim 156$ meV. The
mass range of the effective Majorana neutrino comes from the nuclear matrix
elements calculated by several models. The lower half-life limits by other
groups are longer than the order of $10^{25}$ yr. One needs to develop an
extremely low background radiation detector to search for $0\nu\beta\beta$
events.

##  2 PICOLON Project

PICOLON (Pure Inorganic Crystal Observatory for LOw-energy Neutr(al)ino) is a
research project to search for WIMPs, and double beta decays, especially
$0\nu\beta\beta$ using extremely pure inorganic crystals. We are currently
developing detectors for the WIMP search using NaI(Tl)[21, 20].

We have established a technique to remove radioactive impurities in NaI(Tl)
crystals. Table 1 shows our recent results and comparisons with other groups
applying NaI(Tl) crystal. We calculated the U-series contamination from 226Ra
since this series is divided into five sub-series by long-lived progeny. The
226Ra produces many isotopes, which emit many beta-rays and gamma-rays. The
beta- and gamma-rays from them are supposed to be the primary origin of the
present background.

Table 1: The concentration of RIs in present NaI(Tl) crystals by the groups using NaI(Tl) [22, 23, 19].  Conc. [$\mu$Bq/kg] | DAMA/LIBRA | COSINE | ANAIS | PICOLON
---|---|---|---|---
40K | $<600$ | $515\sim 1900$ | $540\sim 1200$ | $<480$
232Th | $2\sim 31$ | $2.5\sim 35$ | $0.4\sim 4$ | $4.6\pm 1.2$
226Ra | $8.7\sim 124$ | $11\sim 451$ | $2.7\sim 10$ | $8.7\pm 1.5$
210Pb | $5\sim 30$ | $50\sim 3800$ | $740\sim 3150$ | $28\pm 5$

Low background measurements are performed at the Kamioka Underground
Observatory of Tohoku University in Hida City, Gifu Prefecture, Japan. We
installed the latest Ingot # 94 in July 2021 inside a 10 cm thick lead and 5
cm thick oxygen-free copper shield. We flushed pure nitrogen to purge 222Rn in
the shield. NaI(Tl) immediately after the installation has much background due
to cosmogenic isotopes 125I ($T_{1/2}=59.408$ d) and 126I ($T_{1/2}=13.11$ d).
Since background due to these isotopes can be reduced by cooling for several
months after installation, we started low-background measurement in January
2021.

Figure 3: Energy spectrum in low energy region taken by the newest ingot # 94.

The energy spectrum in the low energy region is shown in Figure 3. The
background level between 5 keVee and 10 keVee was as low as 1.4 day-1kg-1keV-1
. Unfortunately, the noise level below 5 keVee was unstable. The noise events
were mainly due to Čerenkov radiation. The second noise origin was a
fluctuation in the baseline of the signal due to insufficient ground
connection. We are trying to remove these noise events by machine-learning.

##  3 CANDLES Project

Figure 4: The side view (a) and top view (b) of the CANDLES-III detector. See
details of the detector setup [24]. Figure 5: The energy spectrum taken for
the live time of 130.4 days [24].

CANDLES (CAlcium fluoride for studies of Neutrino and Dark matters by Low
Energy Spectrometer) aims to search for $0\nu\beta\beta$ of 48Ca by pure CaF2
crystal [24]. We used 96 modules of pure CaF2 crystals, where pure CaF2 means
an undoped CaF2 crystal. The dimension of each CaF2 crystal is 10 cm cubic,
and the total mass of the 48Ca is 350 g. The schematic drawing of the CANDLES-
III detector is shown in Figure 4.

The experiment to search for $0\nu\beta\beta$ event was performed for the live
time of 130.4 days. We selected 21 pure CaF2 modules from 96 ones. The upper
limit on the half-life of $0\nu\beta\beta$ was $5.6\times 10^{22}$ yr at 90 %
C.L. The energy spectrum around the region of interest is shown in Figure 5.
We selected the best 21 crystals from 96 ones because of the higher
contamination of radioactive isotopes.

The background energy spectrum was well reconstructed by known origins. The
background events consist of external gamma-rays from (n,$\gamma$) reaction in
the surrounding rock [25] and internal contamination (green line). The severe
internal contamination is observed from 208Tl and 212Bi-212Po events[24]
around $Q$-value (blue line). The $2\nu\beta\beta$ events also an origin of
background because of the poor enegy resolution (magenta line). The red line
gives the total amount of the known background events. We must reduce these
contaminants lower than 10 $\mu$Bq/kg for further investigation of
$0\nu\beta\beta$.

##  4 Discussion

###  4.1 Prospect of purification

We discuss the prospect of further purification of inorganic crystals. We have
successfully purified NaI(Tl) crystal because the NaI is water-soluble and
thus easy to process. Water-insoluble inorganic crystals, such as CaF2, have
been applied in many double-beta decay experiments. However, we must apply a
different purification method to these crystals. Previously, we tried to wash
CaF2 powder with nitric acid and remove adhered impurities, but this method
cannot clean impurities inside the CaF2 powder.

We started the purification work of CaF2 as a water-insoluble inorganic
crystal. Our target is less than 10 $\mu$Bq/kg for the uranium and thorium
series. In addition, we plan to enrich the target isotope 48Ca in CaF2, so the
purification method wastes as little calcium as possible during the
purification. After considering these conditions, we chose the method of
synthesizing CaF2 after purifying the raw materials of CaF2.

The raw materials of CaF2 are CaCl2 and HF. The HF has achieved an ultra-high
purity of 12 N because it is used in the wet etching process of silicon
wafers. Therefore, we purified CaCl2 by checking its purity.

Before purification, we measured impurities in CaCl2 by ICP-MS. We bought CaF2
powder from three providers, A, B, and C. The concentration of 238U and 232Th
were measured by Agilent 8900 at Tsukuba University. The concentration of 238U
was 166 ppt, 18 ppt, and 30 ppt in samples A, B, and C, respectively.
Unfortunately, the concentration of 232Th was not measured successfully
because of the trouble in ICP-MS.

We tested the purification effectiveness by using a resin that catches uranium
ions. After repeating the purification process three times, the concentration
of 238U was almost equal to that of the blank sample, confirming that the
purification process was successful. In the future, we will produce a high-
purity, large (5 cm cubic) CaF2 crystal and measure its impurities.

###  4.2 Dark matter search by huge inorganic crystal

Current NaI(Tl) detectors consist of small modules whose mass is around 10 kg.
When an array of small modules is used, the background caused by components
other than crystals, such as housings and light guides, increases and cannot
be removed by event screening. XENON-nT is a bulk detector and can remove
external background events by measuring the location of events in the
detector.

It is possible to select events in the center of the detector and remove the
external background by making large inorganic crystals of 50 cm or larger.
Utilizing multiple PMTs gives the precise position of an event occurs. The
ratio of the detected number of scintillation photons enables the removal of
the background events from surrounding materials.

Figure 6: Blue: Before external veto. Red: After external veto.

We estimated the reduction efficiency by installing an outer detector for an
external veto by Monte Carlo simulation (Geant4). The difference in the energy
spectrum before and after the veto is shown in Figure 6. We constructed a
complex NaI(Tl) detector in this simulation. The inner crystal was columnar
shaped with 50 cm diameter and 50 cm length. The inner crystal was covered
with 5 cm thick NaI(Tl) crystal. A PMT was placed on the surface of the outer
NaI(Tl) crystal and radioactive isotopes, 214Pb, 214Bi, and 208Tl were
contained at the divider circuit of the PMT.

We removed the events in that the gamma rays interacted with inner and outer
crystals. The blue histogram is the original energy spectrum, and the red is
the remaining events. A reduction of almost two orders of magnitude was
derived. We are optimizing the detector design and estimating the expected
sensitivity.

##  5 Advantages and disadvantages of inorganic crystals

##  6 PICOLON project

##  7 Discussion

##  8 Acknowledgment

We acknowledge the support of the Kamioka Mining and Smelting Company. This
work was supported by JSPS KAKENHI Grant No. 26104008, 19H00688, 20H05246, and
Discretionary expense of the president of Tokushima University. This work was
also supported by the World Premier International Research Center Initiative
(WPI Initiative).

## References

  * [1] Y. Sofue, Galaxies, 8, 37 (2020).
  * [2] M. Honma et al., Publ. Astron. Soc. Japan, 64, 136 (2012).
  * [3] M. Honma, T. Nagayama, and N. Sakai, Publ. Astron. Soc. Japan, 67, 70 (2015).
  * [4] K. Freese, J. Frieman, and A. Gould, Phys. Rev. D, 37, 3388–3405 (1988).
  * [5] A K. Drukier, K Freese, and D N. Spergel, Phys. Rev. D, 33, 3495–3508 (1986).
  * [6] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett., 38, 1440–1443 (1977).
  * [7] Steven Weinberg, Phys. Rev. Lett., 40, 223–226 (1978).
  * [8] Agnese et al., Phys. Rev. Lett., 112, 241302 (2014).
  * [9] D Cline et al., Astroparticle Physics, 12, 373–377 (2000).
  * [10] E. Aprile and Others, Phys. Rev. Lett., 129, 161805 (2022).
  * [11] H Ejiri, J. Phys. Soc. Japan, 74(8), 2101–2127 (2005).
  * [12] B Jones,  The physics of neutrinoless double beta decay: A beginners guide,  In Proceedings of Theoretical Advanced Study Institute 2020 ”The Obscure Universe: Neutrinos and Other Dark Matters” - TASI2020 — PoS(TASI2020). Sissa Medialab (2021).
  * [13] S. Abe et al., Phys. Rev. Lett., 130, 051801 (2023).
  * [14] D. Abrams et al., Phys. Rev. D, 66, 122003 (2002).
  * [15] R. Agnese et al., Phys. Rev. Lett., 116, 071301 (2016).
  * [16] E. Aprile et al., Eur. Phys. J.C, 78 (2018).
  * [17] E. Aprile et al., Phys. Rev. Lett., 129, 161805 (2022).
  * [18] E. Aprile et al., The European Physical Journal C, 82 (2022).
  * [19] B. J. Park et al., Eur. Phys. J. C, 80, 814 (2020).
  * [20] K. Fushimi et al., Prog. Theor. Exp. Phys., ptab020 (2021).
  * [21] K. Fushimi et al., JPS Conf. Proc., 11, 020003 (2016).
  * [22] R. Bernabei et al., Nucl. Instru. and Meth. in Phys. Res. A, 592, 297–315 (2008).
  * [23] J. Amaré et al., The European Physical Journal C, 79 (2019).
  * [24] S. Ajimura et al., Phys. Rev. D, 103, 092008 (2021).
  * [25] K. Nakajima et al., Astroparticle Physics, 100, 54–60 (2018).


# Protein Structure Prediction Using Basin-Hopping

Michael C. Prentiss,
Center for Theoretical Biological Physics, University of California at San
Diego,
La Jolla, CA 92093, USA
David J. Wales,
University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW,
United Kingdom
and Peter G. Wolynes
Center for Theoretical Biological Physics, University of California at San
Diego,
La Jolla, CA 92093, USA

###### Abstract

Associative memory Hamiltonian structure prediction potentials are not overly
rugged, thereby suggesting their landscapes are like those of actual proteins.
In the present contribution we show how basin-hopping global optimization can
identify low-lying minima for the corresponding mildly frustrated energy
landscapes. For small systems the basin-hopping algorithm succeeds in locating
both lower minima and conformations closer to the experimental structure than
does molecular dynamics with simulated annealing. For large systems the
efﬁciency of basin-hopping decreases for our initial implementation, where the
steps consist of random perturbations to the Cartesian coordinates. We
implemented umbrella sampling using basin-hopping to further conﬁrm when the
global minima are reached. We have also improved the energy surface by
employing bioinformatic techniques for reducing the roughness or variance of
the energy surface. Finally, the basin-hopping calculations have guided
improvements in the excluded volume of the Hamiltonian, producing better
structures. These results suggest a novel and transferable optimization scheme
for future energy function development.

##  1 Introduction

The complexity of the physical interactions that guides the folding of
biomolecules presents a significant challenge for atomistic modeling. Many
current protein models use a coarse grained approach to remove degrees of
freedom, such as non-polar hydrogens, which increases the feasible time step
in molecular dynamics simulations [1, 2]. For a more dramatic improvement of
the computational efficiency solvent degrees of freedom can be reduced [3]. In
this case more severe approximations can prevent the model from reproducing
experimental results. Another option is to reduce the number of degrees of
freedom of the solute. The associative memory Hamiltonian (AMH) [4, 5, 6], is
a coarse-grained molecular mechanics potential inspired by physical models of
the protein folding process, but flexibly incorporates bioinformatic data to
predict protein structure. The AMH is first optimised using the minimal
frustration principle in terms of the $T_{f}/T_{g}$ ratio, which estimates the
separation in energy relative to the variance for the misfolded ensemble.
Along with using the energy of the native structure to estimate $T_{f}$, a
random energy model[7] estimate of the glass transition temperature, $T_{g}$,
is used based on a set of decoy structures. $T_{g}$ represents a
characteristic temperature scale at which kinetic trapping in misfolded states
dominates the dynamics. An improved potential is next obtained that uses
better estimates of the $T_{f}/T_{g}$ ratio obtained by maximizing the
normalized difference between the native state and a sampled set of misfolded
decoys which are self-consistently obtained from the potential itself. The
potential so obtained is transferable for the prediction of structures outside
of the training set. The ratio $T_{f}/T_{g}$ has provided a powerful metric
for the optimisation of both this bioinformatically informed energy function
[8, 9], as well as other types of energy functions incorporating only physical
information [10, 11, 12].

While the optimisation [13] of parameters using a training set of evolved
natural proteins smooths the energy landscape from what it would be for a
random hetero-polymer, the common problem of multiple competing minima
persists even when using a reasonably accurate structure prediction potential,
such as this one. Simulated annealing with molecular dynamics has previously
been used to search the rugged landscapes of optimised structure prediction
potentials [14]. While free energy profiles indicate that better structures
actually are present at low temperatures, the slow kinetics of a glass-like
transition during annealing has prevented these minima from being reached
[15]. To quantitatively investigate the origin of the sampling difficulties it
is desirable to use different search strategies.

Here we implement the basin-hopping global optimisation algorithm [16, 17,
18], which has proved capable of overcoming large energetic barriers in a wide
range of different systems. Basin-hopping is an algorithm where a structural
perturbation is followed by energy minimisation. This procedure effectively
transforms the potential energy surface, by removing high barriers, as shown
in Fig. 1. Moves between local minima are accepted or rejected based upon a
Monte Carlo criterion. Avoiding barriers by employing a numerical minimisation
step not only facilitates movement between local minima, but also broadens
their occupation probability distributions, which overlap over a wider
temperature range, thereby increasing the probability of interconversion [19].
Furthermore, it does not alter the nature of the local minima since the
Hamiltonian itself is not changed, enabling the comparison between molecular
dynamics and basin-hopping generated minima. This method has previously been
applied to find global minima in atomic and molecular clusters [20, 21],
biopolymers [22, 23], and solids [24]. Since the algorithm only requires
coordinates, energies, and gradients, it can be transferred between different
molecular systems such as binary Lennard-Jones clusters, all-atom
biomolecules, or coarse-grained proteins models as in this study.

##  2 Theory and Computational Details

The AMH energy function used in the present work has previously been optimised
over a set of non-homologous $\alpha$ helical proteins, and consists of a
backbone term, $E_{\rm back}$, and an interaction term, $E_{\rm int}$, which
has an additive form [25, 26]. This model is sometimes termed the AMC model
(associative memory contact) to distinguish it from one that uses nonadditive
water mediated interactions termed the AMW model [27, 14]. Since this model
has been described in detail before [15, 28], we will only summarize its form
here. We employ a version of the coarse-grained model where the twenty letter
amino acid code has been reduced to four, and the number of atoms per residue
is limited to three (Cα, Cβ, and O), except for glycine residues. The units of
energy and temperature were both defined during the parameter optimisation.
The interaction energy $\epsilon$ was defined in terms of the native state
energy excluding backbone contributions, $E^{\rm N}_{\rm amc}$, via

$\epsilon=\frac{\left|E^{\rm N}_{\rm amc}\right|}{4N},$ (1)

where $N$ is the number of residues of the protein being considered.
Temperatures are quoted in terms of the reduced temperature ${T_{amc}}=k_{\rm
B}T/\epsilon$. While $E_{\rm back}$ creates self-avoiding peptide-like
stereochemistry, $E_{\rm int}$ introduces the majority of the attractive
interactions that produce folding. The interactions described by $E_{\rm int}$
depend on the sequence separation $\left|i-j\right|$. The interaction between
residues less than 12 amino acids apart were defined by Eqs. (2).

$E_{\rm local}=-\frac{\epsilon}{a}\sum_{\mu=1}^{N_{mem}}\sum_{j-12\leq i\leq
j-3}\gamma(P_{i},P_{j},P_{i^{\prime}}^{\mu},P_{j^{\prime}}^{\mu},x(\left|i-j\right|))\exp\left[-\frac{(r_{ij}-r_{i^{\prime}j^{\prime}}^{\rm\mu})^{2}}{2\sigma^{2}_{ij}}\right],$
(2)

The index $\mu$ runs over all $N_{mem}$ memory proteins to which the protein
has previously been aligned using a sequence-structure threading algorithm
[29] (i.e. each $i$-$j$ pair in the protein has an $i^{\prime}$-$j^{\prime}$
pair associated with it in every memory protein; if, due to gaps in the
alignment, there is no $i^{\prime}$-$j^{\prime}$ pair associated with $i$-$j$
for a particular memory then this memory protein simply gives no contribution
to the interaction between residues $i$ and $j$). The interaction between Cα
and Cβ atoms is a sum of Gaussian wells centred at the separations
$r_{i^{\prime}j^{\prime}}^{\rm\mu}$ of the corresponding memory atoms. The
widths of the Gaussians are given by $\sigma_{ij}=\left|i-j\right|^{0.15}$ Å.
The scaling factor $a$ is used to satisfy Eq. 1. The weights given to each
well are controlled by
$\gamma(P_{i},P_{j},P_{i^{\prime}}^{\mu},P_{j^{\prime}}^{\mu},x(\left|i-j\right|))$,
which depends on the identities $P_{i^{\prime}}$ and $P_{j^{\prime}}$ of the
residues to which $i$ and $j$ are aligned, as well as the identities $P_{i}$
and $P_{j}$ of $i$ and $j$ themselves. The self-consistent optimisation
calculates the $\gamma$ parameter which originates the cooperative folding in
the model. A three well contact potential [Eq. (3)] is used for residues
separated by more than 12 residues as described by,

$E_{\rm
contact}=-\frac{\epsilon}{a}\sum_{i<j-12}\sum_{k=1}^{3}\gamma(P_{i},P_{j},k)c_{k}(N)U(r_{min}(k),r_{max}(k),r_{ij}).$
(3)

Here, the sequence indices $i$ and $j$ sum over all pairs of Cβ atoms
separated by more than 12 residues. The sum over $k$ is over the three wells
which are approximately square wells between $r_{min}(k)$ and $r_{max}(k)$
defined by,

$U(r_{min}(k),r_{max}(k),r_{ij})=\frac{1}{4}\left\\{\left[1+\tanh\left(7[r_{ij}-r_{min}(k)]/\text{\AA}\right)\right]+\left[1+\tanh\left(7[r_{max}(k)-r_{ij}]/\text{\AA}\right)\right]\right\\}.$
(4)

The parameters ($r_{min}(k),r_{max}(k)$), are (4.5Å, 8.0Å) , (8.0Å, 10.0Å),
and (10.0Å, 15.0Å) for $k=1,2$ and 3 respectively. In order to approximately
account for the variation of the probability distribution of pair distances
with number of residues in the protein ($N$) a factor $c_{k}(N)$ has been
included in $E_{\rm long}$. It is given by $c_{1}=1.0$,
$c_{2}=1.0/(0.0065N+0.87)$ and $c_{3}=1.0/(0.042N+0.13)$. The individual wells
are also weighted by $\gamma$ parameters which depend on the identities of the
amino acids involved, using the 4-letter code defined above. In contrast to
the interactions between residues closer in sequence, this part of the
potential does not depend on the database structures that define local-in-
sequence interactions.

To pinpoint the effects of frustration or favorable non-native contacts always
present in any coarse gained protein model, we simulated a perfectly smooth
energy function often called a Gō model [30]. Gō models are an essential tool
for understanding protein folding kinetics [31, 32]. While having the same
backbone terms [33], in this single structure based Hamiltonian [Eq. (5)], all
of the interactions, $E_{\rm int}$ are defined by Gaussians whose minima are
located at the pair distribution found in the experimental structure:

$E_{\rm G{\bar{o}}}^{\rm AM}=-\frac{\epsilon}{a_{\rm G{\bar{o}}}}\sum_{i\leq
j-3}\gamma_{\rm
G{\bar{o}}}(x(\left|i-j\right|))\exp\left[-\frac{(r_{ij}-r_{ij}^{\rm
N})^{2}}{2\sigma^{2}_{ij}}\right].$ (5)

The global minima of such an energy function should be the input structure.

Many have used additional constraining potentials to characterise unsampled
regions of coordinate space while using molecular dynamics [34, 15]. To
characterize the landscape sampled with basin-hopping, we also used a
structure constraining potential to identify ensembles with fixed but varying
fractions of native structure. Using such a potential allows us to access
interesting configurations that are unlikely to be thermally sampled. The
constraining potential also called umbrella potentials are centered on
different values of an order parameter to sample along the collective
coordinates. One of the collective coordinates is $Q$, an order parameter that
measures the sequence-dependent structural similarity of two conformations by
computing the normalized summation of C-alpha pairwise contact differences, as
defined in Eq. (6):[15]

$Q=\frac{2}{(N-1)(N-2)}\sum_{i<j-1}\exp\left[-\frac{(r_{ij}-r_{ij}^{\rm
N})^{2}}{\sigma_{ij}^{2}}\right].$ (6)

The resulting order parameter ranges from zero, where there is no similarity
between structures, to one, which represents an exact overlap. The form of the
potential is $E(Q)=2500\epsilon(Q-Q_{i})^{4}$, where $Q_{i}$ may be varied in
order to sample different regions of the chosen order parameter. As in
equilibrium sampling, simulations were initiated at the native state and the
$Q_{i}$ parameter was reduced throughout the sampling.

We have also studied the potential energy landscape when multiple surfaces are
superimposed upon each other by the use of multiple homologous target
proteins. This manipulation of the energy landscape has been shown to further
reduce local energetic frustration that arises from random mutations in the
sequence away from the consensus optimal sequence for a given structure. By
reducing the number of non-native traps, this averaging often improves the
quality of structure prediction results [35, 36, 37, 38]. As seen in Eq. (7),
the form and the parameters of the energy function are maintained from Eqs.
(2) and (3), but the normalized summation is taken over a set of homologous
sequences:

$E_{\rm
AM}=-\frac{1}{N_{seq}}\sum^{seq}_{k=1}\sum^{N}_{i<j}E_{int}(P^{k}_{i},P^{k}_{j}).$
(7)

Since proteins are not random heteropolymers, the differences in the energy
function for homologous proteins are randomly distributed, therefore the mean
over multiple energy functions should have less energetic variation than the
original function. Indeed, performing this summation is a way of incorporating
to the optimisation of the $T_{f}/T_{g}$ criterion into any energy function.
The target sequences of the homologues can be identified using PSI-Blast with
default parameters [39, 40]. Some classes of proteins have a large number of
sequence homologues, therefore performing a multiple sequence alignment can be
impractical. Removing redundant sequences from within the set of identified
homologues also removes biases that can be introduced where there are few
homologues available. This is done by preventing sequences in the collected
sequences from having greater than 90% sequence identity. The remaining
sequences are aligned in a multiple sequence alignment [41]. Gaps within the
sequence alignment can be addressed within the AMH energy function in a
variety of ways. In the present work, gaps in the target sequence were
removed, while gaps within homologues were completed with residues from the
target protein. While this procedure may introduce small biases toward the
target sequence, it is preferable to ignoring the interactions altogether.

Finally, we made several ad hoc changes to the backbone potential, $E_{\rm
back}$. Eliminating some compromises necessary for rapid molecular dynamics
simulations allowed the AMH potential to be adapted to basin-hopping. Another
goal was to prevent the over-collapse of the proteins by altering the excluded
volume energy term, which should reduce the number of states available during
minimisation. The terms shown in Eq. (8) are used to reproduce the peptide-
like conformations in the original molecular dynamics energy function:

$E_{\rm back}=E_{\rm ev}+E_{\rm harm}+E_{\rm chain}+E_{\rm chi}+E_{\rm Rama}.$
(8)

$E_{\rm ev}$ maintains a sequence specific excluded volume constraint between
the Cα-Cα, Cβ-Cβ, O-O, and Cα-Cβ atoms that are separated by less than $r_{\rm
ev}$. Previously [26], we have seen that modifying $E_{\rm back}$ can produce
a less frustrated energy surface when using thermal equilibrium sampling, but
slow dynamics was often found to result since the local barrier heights became
too large. The ability of basin-hopping to overcome such large, but local
barriers allows us therefore to consider a potential whose dynamics would
otherwise be too slow for molecular dynamics. In the final part of the paper,
we altered the excluded volume term, as shown in Eq. (9) to prevent over-
collapse:

$E_{\rm ev}=\epsilon\lambda_{\rm RV}^{\rm C}\sum_{x,y}\sum_{i<j}\theta(r_{\rm
ev}^{\rm C}(j-i)-r_{{\rm C}^{x}_{i}{\rm C}^{y}_{j}})(r_{\rm ev}^{\rm
C}(j-i)-r_{{\rm C}^{x}_{i}{\rm C}^{y}_{j}})^{2}+\epsilon\lambda_{\rm ev}^{\rm
O}\sum_{i<j}\theta(r_{\rm ev}^{\rm O}-r_{{\rm O}_{i}{\rm O}_{j}})(r_{\rm
ev}^{\rm O}-r_{{\rm O}_{i}{\rm O}_{j}})^{2},$ (9)

by changing the default molecular dynamics parameters, $\lambda_{\rm EV}^{\rm
C}=20$, $\lambda_{\rm EV}^{\rm O}=20$, $r_{\rm ev}^{\rm C}(j-i<5)=3.85\,$Å,
$r_{\rm ev}^{\rm C}(j-i\geq 5)=4.5\,$Å, and $r_{\rm ev}^{\rm O}=3.5\,$Å, to
$\lambda_{\rm EV}^{\rm C}=250$, $\lambda_{\rm EV}^{\rm O}=250$, $r_{\rm
ev}^{\rm C}(j-i<5)=3.85\,$Å, $r_{\rm ev}^{\rm C}(j-i\geq 5)=3.85\,$Å, and
$r_{\rm ev}^{\rm O}=3.85\,$Å. The force constant are over an order of
magnitude larger than those used in the molecular dynamics, and the radii of
the Cα, Cβ, and O atoms are also 10% larger than previously used values. This
increase in excluded volume slows the onset of chain collapse, but improves
steric interactions. The other change to the backbone potential involves terms
which maintain chain connectivity. In molecular dynamics annealing, covalent
bonds are preserved using the SHAKE algorithm [42], which permits an increase
of the molecular dynamics time step. For basin-hopping in all parts of this
paper, we removed the SHAKE method and replaced it with a harmonic potential,
$E_{\rm harm}$, between the Cα-Cα, Cα-Cβ, and Cα-O atoms. This replacement
permits the location of local minima without requiring an internal coordinate
transformation, and avoids discontinuous gradients. When minimised, the
additional harmonic terms typically contribute only only .015 $k_{B}T$ per
bond. The remaining terms of the original backbone potential are maintained.
Depending on the sidechain, the neighbouring residues in sequence sterically
limit the variety of positions the backbone atoms can occupy, as evidenced in
a Ramachandran plot [43]. This distribution of coordinates is reinforced by a
potential, $E_{\rm Rama}$, with artificially low barriers to encourage rapid
local movements. The planarity of the peptide bond is ensured by a harmonic
potential, $E_{\rm chain}$. The chirality of the Cα centres is maintained
using the scalar triple product of neighbouring unit vectors of carbon and
nitrogen bonds, $E_{\rm chi}$.

For basin-hopping simulations, whose algorithm is outlined in Fig. 2 the most
important sampling parameters are the temperature used in the accept/reject
steps for local minima $T_{bh}$, and the maximum step size for perturbations
of the Cartesian coordinates $d$. A higher temperature allows transitions to
an increased energy minima to be accepted, and also creates a larger the
number of iterations typically required to minimise the greater perturbed
configurations. Too high a temperature leads to insufficient exploration of
low-energy regions. The temperature ($T_{bh}$) for these simulations was 10
$T_{amc}$. Lower temperatures resulted in slower escape rates from low energy
traps, while higher temperatures prevented adequate exploration of low energy
regions. The step size needs to be large enough to move the configuration into
the basin of attraction of one local minimum to a neighbouring one, but not be
so large that the new minimum is unrelated to the previous state. Every
Cartesian coordinate was displaced up to a maximum step size($d$) of
$0.75\,$Å, the optimum value determined from preliminary tests. Each run
consisted of 2500 basin-hopping steps saving structures every 5 basin-hopping
steps. The convergence condition ($\delta E_{min}$) on the root-mean-square
(RMS) gradient for each minimisation was set to $10^{-3}\,\epsilon$, and the 5
lowest-lying minima from each run were subsequently converged more tightly
($\delta E_{final}$) to an RMS gradient of $10^{-5}\epsilon$. It is important
to note that basin-hopping does not provide equilibrium thermodynamic
sampling. In structure prediction there, however, is no rigorous need for the
search to obey detailed balance since the global energy minimum is the primary
interest. Basin-hopping provides a means for the optimal global search of the
energy landscape, however other methods must be used when calculating free
energy and entropy.

In previous structure prediction studies with the AMH, low energy structures
were identified using off-lattice Langevin dynamics with simulated annealing,
employing a linear annealing schedule of 10000 steps from a temperature of 2.0
to 0.0, starting from a random configuration [5]. The number and length of
simulations needed in both strategies were determined by the number of
uncorrelated structures encountered. The current basin-hopping method with the
AMC energy function encounters roughly one deep trap per run. In order to
sample 100 independent structures in molecular dynamics 20 separate runs were
needed, because simulated annealing samples about five independent states
before the glass transition temperature is encountered, as measured by the
rapid decay of structural correlations. We compared several $\alpha$ helical
proteins, both from within and outside the training set of the AMH energy
function.

##  3 Results and Discussion

We performed initial calculations with a Gō potential for the 434 repressor
(protein data bank (PDB [44]) ID 1r69). In Fig. 3 we show this model
accurately represents the native basin. Steps where the energy increases are
allowed by the sampling method and are not examples of frustration. Studies on
the Gō model provide a useful benchmark for comparing the computer time
required for the different global optimisation strategies. Using the sampling
parameters used in this report, we compared the time for initial collapse
between the molecular dynamics and basin-hopping runs. The initial collapse
required about 7 minutes for the annealing runs and 31 minutes for basin-
hopping on a desktop computer. However, these values do not reflect the actual
performance of the two approaches in locating global minima, which will depend
upon the move sets, step size, temperature, and convergence criteria.

While using the AMC structure prediction Hamiltonian, we found that basin-
hopping was often able to locate lower energy structures and also identified
minima that have greater structural overlap with the native state than
annealing. These results are produced for structure predictions for proteins
both inside and outside the training set, as demonstrated in Table 1. The
first three proteins (PDB ID 1r69, 3icb, 256b) in Table 1 are in the training
set of the Hamiltonian [25], while the other three are not, and can therefore
be considered as predictions from the algorithm. The minima located with
basin-hopping show an increase in structural overlap with the native state
[Eq. (6)] when compared to the Langevin dynamics approach. $Q$ scores of 0.4
for single domain proteins generally correspond to a low resolution root mean
square deviation (RMSD) of around 5 Åor better. $Q$ scores of 0.5 and higher
have still more accurate tertiary packing and are of comparable quality to the
experimentally derived models. The high quality structures obtained suggest
the form of the backbone terms is appropriate, since the physically correct
stereochemistry is reproduced. Lower energy structures are sampled by basin-
hopping for the non-training set proteins, but the structural overlap
improvement found in these deeper minima was smaller. Larger proteins pose a
greater challenge for basin-hopping with this Hamiltonian due to the random
steps in Cartesian coordinates. Dihedral coordinate move steps would probably
be more efficient, and will be considered in future work.

The distribution of minima encountered from multiple simulations for both
search methods is shown in Fig. 4 where a greater density of high quality
structures is obtained by the basin-hopping algorithm. The potential energy
surface still includes, therefore, significant residual frustration in the
near-native basin in the form of low-lying minima separated by relatively high
barriers. Without the parameter optimisation to reduce frustration, folding
would exhibit more pronounced glassy characteristics. Most of the cooperative
folding occurs during collapse until $Q$ values of around 0.4 are reached.
While the structures from simulated annealing are accurate enough for
functional determination, we see basin-hopping can better overcome barriers
that are created after collapse. The density of the high quality structures is
also important for post-simulation k-means clustering analysis [45]. Another
way of representing the data of a set of independent basin-hopping simulations
is by selecting the lowest energy structures from each simulation of the 434
repressor (PDB ID 1r69) and HDEA (PDB ID 1bg8) proteins and ordering them with
respect to their structural overlap. As shown in Fig. 5, the protein in the
training set (434 repressor) produces better results than the non-training
protein, as expected.

We have decomposed the different energy terms in the Hamiltonian in Table 2,
to examine which interactions are most effectively minimised. The AMC
potential has three different distance classes in terms of sequence
separation, and these are defined as short ($\left|i-j\right|<5$), medium
($5\leq\left|i-j\right|\leq 12$), and long ($\left|i-j\right|>12$). Most
importantly, the long-range AMH interactions are successfully minimised in the
basin-hopping runs, due to the ability of basin-hopping to overcome large
energetic barriers. This term will govern the quality of structures sampled
using an approximately smooth energy landscape. The other terms that define
secondary structure formation are not as well minimised. This result is due to
the disruption of helices by the random Cartesian perturbation move steps.
These move steps benefit favorable steric packing and therefore do well at
minimising the excluded volume energy term of the Hamiltonian. A combined
minimisation approach might be more efficient, where larger dihedral steps
could be made early during minimisation to sample a wider number of
structures, followed by random Cartesian steps to optimise the steric
interactions.

While we sampled high quality structures, we would like to confirm that we
have completely sampled the global minima of the energy surface. To access
these unsampled states we used umbrella potentials. When constraining a set of
simulations to different values of $Q$, we have obtained energy minima for
cytochrome c, roughly $15\,\epsilon$ deeper than those from unconstrained
minimisations starting with a randomized structure, as shown in Fig. 6. For
the 434 repressor the minima obtained from randomized states and those found
with the $Q$ constraints applied differ by only a few $k_{B}T$. This shows
basin-hopping does indeed act as a global optimisation method, by accurately
identifying the global energy minimum from multiple independent unconstrained
simulation. This behavior is predictable from the choices that governed the
design of the Hamiltonian. Low energy barriers between structures are
desirable during a molecular dynamics simulation because they accelerate the
dynamics. However, for basin-hopping these low barriers encourage tertiary
contact formation before secondary structure units condense for sequences
greater than 110 amino acids.

### Superposition of Multiple Energy Landscapes

Constructing a Hamiltonian by calculating the arithmetic average of the
potential over a set of homologous sequences increased the quality of
predictions in both equilibrium and annealing simulations. We have found this
approach also improved the performance in basin-hopping simulations. For two
different proteins, 100 independent basin-hopping runs were performed with
both the standard and sequence-averaged Hamiltonians. By the superposition of
multiple energy landscapes we saw a reduction in the number of competing low
energy traps around $Q$ values of $0.3$ for both the 434 repressor and
uteroglobin (PDB ID 1UTG), as shown in Fig. 7. Improvement of structure
prediction Hamiltonians can be statistically described by the average energy
gap between the native basin and a set of unfolded structures, and by the
roughness of the energy surface, which corresponds to the variance of the
energy. The sequence based energy function summations limited the energetic
variance of the sampled landscapes, thereby reducing the glass transition
temperature. This improvement, even at the low temperatures sampled in basin-
hopping, is predicted from theory, but difficult to observe in conventional
equilibrium simulations due to the emergent glassy dynamics, which slows the
kinetics. The energy gap improvement was smaller than the reduction of the
energetic variation of the Hamiltonian. In terms of the goal of maximizing the
ratio of $T_{f}/T_{g}$, this increase came primarily from to reducing the
glass transition temperature $T_{g}$. In the low energy region we saw fewer
competing states, and an increased correlation between ${E}$ and ${Q}$ for the
sequence-averaged Hamiltonian compared to the original Hamiltonian. For the
434 repressor the lowest energy structure had the highest $Q$ value
encountered.

### Characterisation of Polymer Collapse

When we annealed the Hamiltonian using molecular dynamics we observed some
over-collapse of the polypeptide chain, producing a smaller radius of gyration
than the experimental structure. In basin-hopping runs we also found
structures exhibiting a larger number of contacts than the experimental
structure, as show in Fig. 8, where a contact is defined as a Cα-Cα distance
of less than 8 Å. While the low-energy structures may be native-like, these
structures were more compact than those observed experimentally. To
investigate this behavior, we examined the backbone and interaction terms of
the Hamiltonian separately using the Gō Hamiltonian in Eq. (5). Somewhat
surprisingly, the Gō model also produces over-collapse, as shown in Fig. 9.
Hence the interaction parameters of the structure prediction Hamiltonian were
not responsible for all of the over-collapse. These minimal model-dependent
frustrations were only eliminated in the final stages of minimisation. The
most effective technique for reducing over-collapse was to increase the force
constant and the atomic radius in the excluded volume terms [Eq. (9)]. The
barrier crossing capabilities of basin-hopping steps produce more over-
collapse than do the annealing minimisations without these parameter changes.
The glass-like transition seen in simulated annealing prevents further
collapse in molecular dynamics, as the rearrangement rates slow down
exponentially with temperature. The improved parameter set of Fig. 10 shows
more native-like collapse, but the lowest energy structures had $Q$ values of
0.36 and the best $Q$ value was 0.45, which are worse than basin-hopping
simulations with the original parameters.

##  4 Conclusion

In this report we have demonstrated that minima with lower energy and higher
quality structures can often be located for the AMH potential using basin-
hopping global optimisation compared to annealing. Encouragingly, the long
range in sequence energy contributions are better minimised than with
simulated annealing. Umbrella sampling using basin-hopping can also show when
the global minima are reached for a selected order parameter. Previous
techniques for reducing the energetic variance of the energy surface in
simulated annealing are also applicable to basin-hopping. Using basin-hopping
also permits improvements in certain backbone terms of the Hamiltonian. These
changes would make the kinetics too slow in molecular dynamics annealing runs,
but larger barriers can be easily crossed using basin-hopping.

These results highlight future optimisation strategies where the deep non-
native traps found by basin-hopping could be used as decoys for further
parameter refinement, rather than the higher-lying minima obtained by
quenching with simulated annealing. This re-optimisation of the potential
results in getting a better estimate for $T_{f}/T_{g}$ now possible because of
the efficiency of the basin-hopping algorithm at identifying low energy
decoys. Another future direction would be evaluating the equilibrium
properties of low-lying structures identified by basin-hopping also provides a
means to calculate free energy barriers, which would be difficult to
characterise via conventional simulations.

## Acknowledgements

We thank Dr Joanne Carr, and Dr Justin Bois for helpful comments throughout
this research. The efforts of P.G.W. and M.C.P. are supported through the
National Institutes of Health Grant 5RO1GM44557. Computing resources were
supplied by the Center for Theoretical Biological Physics through National
Science Foundation Grants PHY0216576 and PHY0225630. M.C.P. gratefully
acknowledges the support by the International Institute for Complex Adaptive
Matter (ICAM-I2CAM) NSF Grant DMR-0645461.

## References

  * [1] J. C. Phillips et al.,  J. Comp. Chem. 26, 1781 (2005).
  * [2] D. V. D. Spoel et al.,  J. Comp. Chem. 26, 1701 (2005).
  * [3] Y. Levy and J. N. Onuchic,  Annu. Rev. Biochem. 35, 389 (2006).
  * [4] M. S. Friedrichs and P. G. Wolynes,  Science 246, 371 (1989).
  * [5] M. S. Friedrichs and P. G. Wolynes,  Tet. Comp. Meth. 3, 175 (1990).
  * [6] M. S. Friedrichs, R. A. Goldstein, and P. G. Wolynes,  J. Mol. Biol. 222, 1013 (1991).
  * [7] B. Derrida,  Phys. Rev. B. 24, 2613 (1981).
  * [8] R. A. Goldstein, Z. A. Luthey-Schulten, and P. G. Wolynes,  Proc. Natl. Acad. Sci. USA 89, 4918 (1992).
  * [9] P. Barth, J. Schonbrun, and D. Baker,  Proc. Natl. Acad. Sci. USA 104, 15682 (2007).
  * [10] J. Lee et al.,  J. Phys. Chem. B 105, 7291 (2001).
  * [11] Y. Fujitsuka, S. Takada, Z. A. Luthey-Schulten, and P. G. Wolynes,  Proteins 54, 88 (2004).
  * [12] Y. Fujitsuka, G. Chikenji, and S. Takada,  Proteins 62, 381 (2006).
  * [13] In this paper we note, the word optimisation will be used with two similar but distinct ways. The first use is for the calculation of the best set of parameters to define a minimally frustrated energy function. The second use is to find the global energy minima on a surface that has multiple minima of nearly equal energies.
  * [14] C. Zong, G. Papoian, J. Ulander, and P. Wolynes,  J. Am. Chem. Soc. 128, 5168 (2006).
  * [15] M. P. Eastwood, C. Hardin, Z. Luthey-Schulten, and P. G. Wolynes,  I.B.M. Systems Research 45, 475 (2001).
  * [16] D. Wales and J. Doye,  J. Phys. Chem. A 101, 5111 (1997).
  * [17] Z. Li and H. A. Scheraga,  Proc. Natl. Acad. Sci. USA 84, 6611 (1987).
  * [18] D. J. Wales and H. A. Scheraga,  Science 285, 1368 (1999).
  * [19] J. P. K. Doye and D. J. Wales,  Phys. Rev. Lett. 80, 1357 (1998).
  * [20] D. Wales,  Energy Landscapes,  Cambridge University Press, Cambridge, UK, 2003.
  * [21] D. J. Wales et al.,  The Cambridge Cluster Database, URL http://www-wales.ch.cam.ac.uk/CCD.html (2001).
  * [22] J. M. Carr and D. J. Wales,  J. Chem. Phys. 123, 234901 (2005).
  * [23] A. Verma, A. Schug, K. H. Lee, and W. Wenzel,  J. Chem. Phys. 124, 044515 (2006).
  * [24] T. F. Middleton, J. Hernández-Rojas, P. N. Mortenson, and D. J. Wales,  Phys. Rev. B 64, 184201 (2001).
  * [25] C. Hardin, M. Eastwood, Z. Luthey-Schulten, and P. G. Wolynes,  Proc. Natl. Acad. Sci. USA 97, 14235 (2000).
  * [26] M. Eastwood, C. Hardin, Z. Luthey-Schulten, and P. G. Wolynes,  J. Chem. Phys. 117, 4602 (2002).
  * [27] G. Papoian, J. Ulander, M. Eastwood, Z. Luthey-Schulten, and P. G. Wolynes,  Proc. Nat. Acad. Sci. U.S.A. 101, 3352 (2004).
  * [28] M. C. Prentiss, C. Hardin, M. Eastwood, C. Zong, and P. G. Wolynes,  J. Chem. Ther. Comp. 2:3, 705 (2006).
  * [29] K. K. Koretke, Z. Luthey-Schulten, and P. G. Wolynes,  Protein Sci. 5, 1043 (1996).
  * [30] N. Gō,  Ann. Rev. Biophys. Bioeng. 12, 183 (1983).
  * [31] J. J. Portman, S. Takada, and P. G. Wolynes,  Phys. Rev. Lett. 81, 5237 (1998).
  * [32] N. Koga and S. Takada,  J Mol Biol 313, 171 (2001).
  * [33] M. P. Eastwood and P. G. Wolynes,  J. Chem. Phys. 114, 4702 (2002).
  * [34] X. Kong and C. L. Brooks III,  J. Chem. Phys. 105, 2414 (1996).
  * [35] F. R. Maxfield and H. A. Scheraga,  Biochemistry 18, 697 (1979).
  * [36] C. Keasar, R. Elber, and J. Skolnick,  Folding and Design 2, 247 (1997).
  * [37] R. Bonneau, C. E. M. Strauss, and D. Baker,  Proteins 43, 1 (2001).
  * [38] C. Hardin, M. P. Eastwood, M. C. Prentiss, Z. Luthey-Schulten, and P. G. Wolynes,  Proc. Natl. Acad. Sci. USA 100, 1679 (2003).
  * [39] S. Altschul et al.,  Nucleic Acids Res. 25, 3389 (1997).
  * [40] J. E. Stajich et al.,  Genome Res. 12, 1611 (2002).
  * [41] J. Thompson, D. Higgins, and T. Gibson,  Nucleic Acids Res. 22, 4673 (1994).
  * [42] J. Ryckaert, G. Ciccotti, and H. Berendsen,  J. Comput. Phys. 23, 327 (1977).
  * [43] G. Ramachandran and V. Sasisekharan,  Adv. Protein. Chem. 23, 283 (1968).
  * [44] H. M. Berman et al.,  Nucleic Acids Res. 28, 235 (2000).
  * [45] D. Shortle, K. T. Simons, and D. Baker,  Proc. Natl. Acad. Sci. USA 95, 11158 (1998).

## Figure Captions

  1. 1.

In the basin-hopping approach the original potential energy surface (solid) is
transformed into a set of plateaus (dashed). The local minima are not changed,
but the transition state regions are removed.

  2. 2.

The basin-hopping algorithm is defined by a few parameters that allow its
transfer to different systems.

  3. 3.

Variation of the energy of the current minimum as a function of $Q$ for minima
encountered in the Markov chain during a basin-hopping run using a Gō model.
Steps that increase the energy are sometimes allowed by the Monte Carlo
criterion, which employed a temperature of 10 $k_{b}$/$\epsilon$.

  4. 4.

Energy as a function of $Q$ for local minima of 434 repressor encountered
during 100 independent basin-hopping optimisations (top) and 20 annealing
simulations (bottom).

  5. 5.

The lowest energy structures of the training set protein, 434 repressor (top)
and the blind prediction proteins, HDEA (bottom) identified from 100
independent basin-hopping simulations. Each minimum has values for energy
illustrated by the dots and structural overlap to the native state $Q$
represented by lines. These minima are ordered with respect to their
structural overlap $Q$ with the native state (Index). The data shows
correlations between the energy and $Q$, while the number of high quality
structures is superior for the training protein.

  6. 6.

Energy as a function of $Q$ for the 434 repressor and cytochrome c proteins
obtained in basin-hopping calculations with the structure prediction
Hamiltonian. These runs employed an additional umbrella potential that
constrains the simulation to different values of $Q$. The results for the 434
repressor are similar to the unconstrained basin-hopping results, but the
structures for cytochrome c are $15\epsilon$ lower in energy than those found
in unconstrained basin-hopping runs.

  7. 7.

Energies of local minima obtained using basin-hopping with the original and a
sequence-averaged Hamiltonian for two training proteins. Importantly for both
the top graph (434 repressor) and the bottom graph (uteroglobin) fewer non-
native states are seen with the sequence averaged (red) Hamiltonian when
compared to standard Hamiltonian (black).

  8. 8.

Results of 100 independent basin-hopping runs for the 434 repressor using the
set of backbone parameters that was optimised for molecular dynamics.
Structures were saved every 20 basin-hopping steps. The ratio of contacts to
native state contacts shows that most of the structures are more compact than
the native state.

  9. 9.

A Gō potential simulation for the 434 repressor shows a modest amount of over-
collapse during a basin-hopping simulation, which is resolved as the structure
approaches a $Q$ value of 1.0.

  10. 10.

Results of 100 independent basin-hopping runs for the 434 repressor using the
set of backbone parameters that was optimised for molecular dynamics.
Structures were saved every 20 basin-hopping steps. An altered set of backbone
parameters produces structures that have similar collapse behavior when
compared to the native state.

## Figures

Figure 1:

  * Basin Hopping Algorithm

  * Monte Carlo Step ($n$ steps)

  * random Cartesian move step with maximum distance ($d$) and temperature ($T_{bh}$)

  * Minimisation

  * L-BFGS quasi-Newtonian method for optimization

  * convergence condition ($\delta E_{min}$) is RMS gradient of $10^{-3}\,\epsilon/r$

  * Minimisation with tight convergence (after $n$ steps)

  * convergence condition ($\delta E_{final}$) is RMS gradient of $10^{-5}\,\epsilon/r$

Figure 2:

Figure 3:

Figure 4:

Figure 5:

Figure 6:

Figure 7:

Figure 8:

Figure 9:

Figure 10:

## Tables

Table 1:  Minima located by molecular dynamics/annealing (MD) and basin-
hopping (BH); the first three proteins are in the training set of the
Hamiltonian, while the results for the second three proteins are predictions.
|  | MD |  |  |  | BH |  |  |
---|---|---|---|---|---|---|---|---|---
PDB ID | length | Lowest $E$ | $Q$ | Highest $Q$ | $E$ | Lowest $E$ | $Q$ | Highest $Q$ | $E$
1r69 | 63 | -428.92 | 0.39 | 0.53 | -307.96 | -435.82 | 0.39 | 0.52 | -408.482
3icb | 75 | -536.98 | 0.47 | 0.52 | -390.54 | -546.57 | 0.40 | 0.49 | -518.92
256b | 106 | -735.02 | 0.42 | 0.65 | -707.51 | -737.31 | 0.37 | 0.40 | -716.51
1uzc | 69 | -457.55 | 0.36 | 0.42 | -383.08 | -458.09 | 0.37 | 0.45 | -433.41
1bg8 | 76 | -469.49 | 0.25 | 0.34 | -465.19 | -468.67 | 0.36 | 0.39 | -461.50
1bqv | 110 | -737.91 | 0.21 | 0.27 | -441.92 | -764.20 | 0.23 | 0.27 | -481.22

Table 2:  Contribution of different energy terms in local minima obtained
using molecular dynamics/annealing (MD) and basin-hopping (BH).
PDB | Method | length | Ex Vol | Rama | Short Range | Medium Range | Long Range
---|---|---|---|---|---|---|---
1r69 | MD | 63 | 9.77 |  $-$101.64 |  $-$128.90 |  $-$84.87 |  $-$123.28
1r69 | BH | 63 | 2.65 |  $-$91.06 |  $-$125.04 |  $-$84.80 |  $-$137.57
3icb | MD | 75 | 11.74 |  $-$127.70 |  $-$177.21 |  $-$90.11 |  $-$153.69
3icb | BH | 75 | 4.40 |  $-$115.76 |  $-$178.47 |  $-$83.37 |  $-$173.38
1uzc | MD | 69 | 10.10 |  $-$118.66 |  $-$134.00 |  $-$90.75 |  $-$124.24
1uzc | BH | 69 | 2.22 |  $-$106.20 |  $-$137.95 |  $-$92.40 |  $-$123.77
1bg8 | MD | 76 | 11.68 |  $-$136.39 |  $-$173.45 |  $-$94.40 |  $-$76.94
1bg8 | BH | 76 | 2.72 |  $-$112.13 |  $-$151.95 |  $-$94.23 |  $-$113.09