# k-Nearest Neighbor Adaptive Sampling, a Simple Tool to Efficiently Explore Conformational Space

Evianne Rovers, Anvith Thudi, Jérôme Hénin, Chris Maddison, Matthieu Schapira

► To cite this version:

Evianne Rovers, Anvith Thudi, Jérôme Hénin, Chris Maddison, Matthieu Schapira. k-Nearest Neighbor Adaptive Sampling, a Simple Tool to Efficiently Explore Conformational Space. Journal of Chemical Theory and Computation, 2025, 21 (17), pp.8267-8277. 10.1021/acs.jctc.5c00462. http://dx.doi.org/10.1021/acs.jctc.5c00462

HAL Id: hal-05241482

https://hal.science/hal-05241482v1

Submitted on 4 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
k-Nearest Neighbour Adaptive Sampling (kNN-AS),
a Simple Tool to Efficiently Explore Conformational
Space

Evianne Rovers¹,²,⁴,*, Anvith Thudi³,⁴, Jérôme Hénin⁵, Chris Maddison³,⁴, Matthieu Schapira¹,²,\*

¹ Structural Genomics Consortium, Toronto, M5G 1L7, Canada

² Department of Pharmacology, University of Toronto, Toronto, M5G 1L7, Canada

³ Department of Computer Science, University of Toronto, Toronto, M5G 1L7, Canada

⁴ Vector Institute, University of Toronto, Toronto, M5G 0C6, Canada

⁵ Laboratoire de Biochimie Théorique UPR 9080, Université Paris Cité, CNRS, 75005 Paris,
France.

* matthieu.schapira@utoronto.ca

* evianne.rovers@mail.utoronto.ca

Abstract

Molecular dynamics (MD) simulations are computationally expensive, a limiting factor when
simulating biomolecular systems. Adaptive sampling approaches can accelerate the exploration of
conformational space by running repeated short MD simulations from well-chosen starting points.
Existing approaches to adaptive sampling have been optimized to either guide sampling in a
desired direction or explore well-formed convex spaces. Here, we describe a novel adaptive
sampling algorithm that leverages a k-nearest neighbour (k-NN) graph of the sampled
conformational space to preferentially launch explorations from boundary states. We term this
approach k-Nearest Neighbor Adaptive Sampling (kNN-AS) and show it has state-of-the-art
performance on simple and complex artificial energy functions and generalizes well on a protein
test case. Implementation of kNN-AS is light, simple and suited to complex real-world applications
where the dimension and shape of the energy landscape is unknown.
Introduction

Molecular dynamics (MD) simulations are used to understand the dynamics of
biomolecular systems, which in turn can help to investigate the molecular mechanism of disease-
related mutations¹⁻⁶, uncover cryptic ligand binding sites⁷⁻¹¹, or predict how allosteric compounds
induce conformational shifts¹²,¹³. MD simulations are computationally intensive, and large systems
often fail to reach equilibrium, let alone explore the full space of conformational states. Indeed,
simulating processes such as protein folding, or protein complex formation, occur on timescales
that can require months or even years of GPU time.

Methods have emerged that seek to explore the free energy landscape of a system more
efficiently and are called enhanced sampling methods¹⁴. The methods can be roughly divided into
adaptive sampling and biased approaches¹⁴. Protocols that introduce a bias, like simulated
annealing and metadynamics¹⁵, alter the energy of the entire system, to try and accelerate
exploration of conformational space¹⁴,¹⁵. However, for large molecular systems these methods
become computationally resource intensive¹⁴. Additionally, natural dynamics of the system is
disturbed which could lead to artifacts when calculating the dynamics from these simulations¹⁵
Lastly, heating a multiprotein system can result in disassembly of the complex rather than
conformational searching as seen for PROTAC-induced ternary complexes¹⁶. Adaptive sampling,
in contrast, seeks to accelerate exploration of conformational space under normal conditions¹⁷.
Here, long MD-simulations are divided into repeated short simulations, where starting states are
carefully selected to optimize for exploration. Algorithms that select new starting states seek to
balance exploration and exploitation of the free-energy landscape defined by a feature set, termed
collective variables (CVs). Since no biasing force was applied to any individual trajectory
simulation, adaptive sampling methods can use Markov State Models (MSM) to construct
equilibrium ensembles despite off-equilibrium initializations¹⁸. Nevertheless, the choice of algorithm is key to ensure that the conformational space is searched efficiently.

Where a specific goal in the conformational space is sought, the gradient of a feature / CV can be used to drive exploration. The algorithms FAST¹⁹ and AdaptiveBandit²⁰ both excel at this task. However, these methods cannot be used when the desired direction or space of the exploration are unknown. REAP/MAREAP use reinforcement learning to learn efficient exploration directions in the features/CVs space²¹. It treats features/CVs independently and rewards states that result in CVs that deviate furthest from the population of known states²¹. The MAREAP algorithm was shown to outperform AdaptiveBandit in a protein conformational sampling task²¹. Nonetheless, this method performs best with linear landscapes where there is only one optimal direction of exploration or when the direction switches from one CV to another (for example, a L-shaped landscape). This is similar to PCA-based sampling which would only work if the landscape can be approximated linearly²². Multi-agent TSLC, a different algorithm, performs better with non-linear landscapes²¹. Since the shape of the landscape is typically unknown, which algorithm is most appropriate *a priori* is unclear. Therefore, an ideal algorithm would efficiently explore the conformational space regardless of dimension and shape. The iMapD²³ and OFLOOD²⁴,²⁵ methods are methods that address these limitations by sampling from nonlinear landscapes. iMapD specifically calculates the underlying manifold of the free energy landscape and uses boundary algorithms to identify points on the edge for new simulations²³. The current implementation is limited to free-energy landscape dimension of 3 due to the boundary detection algorithms, but could be extended with algorithms that could detect boundaries in higher dimension. The algorithm described in this paper, for example, could be explored in conjunction with the manifold learning
approach of iMapD. OFLOOD starts from outliers detected in the free energy landscape by clustering the conformation space and determining the global density²⁴,²⁵. To further explore underrepresented regions, simulations are initiated from clusters that are less populated. Another method could be to use a k-Nearest Neighbour (k-NN) graph to identify nodes located at the boundary of the graph by either comparing the average distance or summing up the vectors towards its neighbours.

In this paper, we present a simple algorithm that leverages a k-Nearest Neighbor (k-NN) graph applied to a set of known molecular states which identifies the edges of the conformational space for focused sampling. We term this approach k-Nearest Neighbor Adaptive Sampling (kNN-AS). This approach builds on extensive work showing the power of k-NN graphs for searching complex feature spaces, for example, as employed by Google's PageRank algorithm²⁶. We briefly describe the algorithm and then demonstrate that it achieves superior performance on both simple and complex artificial energy functions, as well as a real-world protein MD simulation case. The code for the algorithm and experiments can be found on github: https://github.com/ERovers/kNN-AS.git

## Algorithm

One goal of adaptive sampling is to select starting states which maximize the number of new states observed during the simulation. The key challenge is finding a local heuristic that globally optimizes this goal. Reinforcement learning (RL) is one approach, but this can introduce modeling biases and/or require many simulations to learn an adequate model. We instead develop a simple heuristic based on our observations of molecular dynamics simulations: we found that the
average magnitude of the sum of vectors from a state to its nearest neighbours is an effective
indicator of how many new states the dynamics will explore (Figure 1). In other words, states on
the edge the free energy landscapes which has states neighbouring on only one side will have a
bigger sum of vectors and will explore more new states (Figure 1). One can also see this score as
approximately picking states with the largest free energy gradient, with the observation that these
coincide with boundary states. In particular, following the first-order taylor approximation of the
gradient (in particular equation 8) used in Carli et al.²⁷, we have that the gradient is approximately
proportional to the sum of differences to the k-nearest neighbours, where we assume k is a constant.
Hence, ordering datapoints by the magnitude of the sum of differences to their neighbours is
approximately the same as ordering them using the magnitude of their free energy gradient.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Points with high overall directional vector (e.g. points on the edge of a point cloud) discover more new states than points with small overall directional vector (e.g. points within point cloud). A) Visualization of how overall directional vector is calculated from a kNN graph. B) Current states identified in MD simulations are plotted in collective variable space and the overall directional vector for each point is calculated using 10 of its closest neighbours. 30 starting points with the largest magnitude of directional vector, around the mean, and the smallest overall directional vectors are selected and the distributions of magnitudes for each selection are shown in</figcaption></figure>
the last plot. Parallel simulations are run with starting points and number of new states are counted where a new state is defined as a point outside of the alphashape obtained from the previous generation point cloud.

So, our adaptive sampling algorithm uses a k-NN approach to construct a graph from a point cloud (states in feature/CV space). The edges that connect a vertex (state) to its k-NNs are summed and new points for the next generation of simulations are then selected based on the magnitude of the sum of vectors (Algorithm 1). Intuitively, this algorithm promotes exploration of states at the boundary of a free energy landscape (Figure 1).

**Algorithm 1.** kNN-AS sampling algorithm to determine states at the boundary of a free energy landscape using sum of vectors as a metric.
<caption>**Algorithm 1:** K-nearest neighbour adaptive sampling (sum of vectors)</caption>

```c
    Data: starting point positions SP, number of neighbours k, number of states to
          select N, number of epochs E, empty vector for summation C, array of saved
          positions POS, array of mapped collective variables based on positions CV,
          P percentage taken for subsampling, CV2 subsample of CV, POS2
          subsample of positions

1  , for e ∈ {1...E} do
2      for n ∈ {1...N} do
3          pos ← MD simulation with SPn;
4          cv ← Project generated pos into collective variable space;
5          POS ← POS ∪ pos;
6          CV ← CV ∪ cv;
7      end
8  I = RandomSubset(values = {1, ..., |CV|}, size = P|CV|);
9  CV2, POS2 ← CV[I], POS[I];
10 distances, indexes ← NearestNeighbours(CV2, k);
11 Scores ← ∅;
12 for cv ∈ CV2 do
13     C ← 0; // all zeros vector with same dimensionality as CV space
14     for i ∈ indexes do
15         C ← C + (cv − CV2i)
16     end
17     scores : scores ∪ ||C||²;
18 end
19 SP ← POS2[argsort(scores)[-N :]]
20 end
```

Another metric to detect points on the boundary would involve measure local density by, for example, calculating the average distance to k-nearest neighbours. The key difference between using the sum of vectors or average distance as metrics for identifying boundary points lies in what they reveal about the landscape. The sum of vectors tends to highlight corner regions, which may indicate potential transition pathways (Figure 2A). In contrast, average distance also captures sparsely sampled areas, which could represent isolated patches within the broader landscape (Figure 2B). Therefore, we also tested an algorithm by generating k-NN graph from a point cloud and measure the distances its k-neighbours and calculate an average of these
distances. The states with the highest average distance to its neighbours are chosen for the next generation of simulations (Algorithm 2).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Conceptual difference between sum of vectors or average distance to neighbours for selecting states on the boundary. A) Sum of vectors identifies directionality in the landscape and thereby potential transition pathways to unexplored areas. B) Average distance to neighbours identifies sparse areas within the landscape and could fill up underexplored areas.</figcaption></figure>

**Algorithm 2.** kNN-AS sampling algorithm to determine states at the boundary of a free energy landscape using average distance to neighbours as a metric.
**Algorithm 2:** K-nearest neighbour adaptive sampling (average distance)

**Data:** starting point positions *SP*, number of neighbours *k*, number of states to select *N*, number of epochs *E*, empty vector for summation *C*, array of saved positions *POS*, array of mapped collective variables based on positions *CV*, *P* percentage taken for subsampling, *CV2* subsample of *CV*, *POS2* subsample of positions

1, **for** *e* ∈ {1...*E*} **do**
2 &nbsp;&nbsp;&nbsp;&nbsp; **for** *n* ∈ {1...*N*} **do**
3 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; *pos* ← MD simulation with *SP*<sub>*n*</sub>;
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; *cv* ← Project generated *pos* into collective variable space;
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; *POS* ← *POS* ∪ *pos*;
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; *CV* ← *CV* ∪ *cv*;
7 **end**
8 &nbsp;&nbsp;&nbsp;&nbsp; *I* = RandomSubset(values = {1, ..., |*CV*|}, size = *P*|*CV*|);
9 &nbsp;&nbsp;&nbsp;&nbsp; *CV2*, *POS2* ← *CV*[*I*], *POS*[*I*];
10 &nbsp;&nbsp;&nbsp;&nbsp; distances, indexes ← NearestNeighbours(*CV2*, *k*);
11 &nbsp;&nbsp;&nbsp;&nbsp; Scores ← ∅;
12 &nbsp;&nbsp;&nbsp;&nbsp; **for** *cv* ∈ *CV2* **do**
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; *C* ← 0
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **for** *i* ∈ indexes **do**
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; *C* ← *C* + ||(*cv* − *CV2*<sub>*i*</sub>)||<sup>2</sup>
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **end**
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; scores : scores ∪ C/k;
18 **end**
19 &nbsp;&nbsp;&nbsp;&nbsp; *SP* ← *POS2*[argsort(scores)[-*N* :]]
20 **end**

## Results

### Sensitivity to hyperparameters

The kNN-AS algorithm has two hyperparameters specific to the algorithm: the size of subsample (P) used for generating the k-NN graph, and number of neighbours (k). When performing a test study on an artificial energy landscape (Figure 3A) using the sum of vectors metric, we found that the size of subsample is the most important hyperparameter, since the extent of exploration increases when the size of the subsample decreases (Figure 3B). The reason for this behaviour could be that subsampling results in a neighbourhood structure of the k-NN graph that
reflects the global structure of a point cloud, as opposed to the local structure induced by the autocorrelations in MD trajectories when using the whole dataset. In particular, during MD simulations, frames are saved every few timesteps. If this interval is too short, the frames can form a nearly continuous line in the conformational landscape, leading k-NN graph to pick neighbours from the same trajectory and causing the graph to only detect local structure. To check for autocorrelation, a global set of angles between the vertices can be calculated which should show a uniform distribution (Figure 4, top right). If a skewed distribution is observed with peaks around 0 and 180 degrees (Figure 4 top left), then k-NN graph only detects the local structure of the MD trajectory itself, namely all points form a line (Figure 4, bottom left). Subsampling addresses this by removing the neighbouring frames, creating a sparser and more random dataset. This results in selecting states on the edges of the explored space (Figure 4, bottom right), which further promotes exploration. Lastly, the state selection step (including generating the kNN graph) can be computationally inexpensive for medium size datasets (<10⁵ datapoints takes <1 second), but for larger datasets we do see an exponential increase in computation time (Figure S1A). Therefore, subsampling also helps with lowering the computational time. The dimensions of a dataset also influence computational time (Figure S1B), where the algorithm takes 1 sec for less than 20 dimensions and increases linearly after. For dimension of 200, the algorithm only takes 3 seconds to perform state selection step. So, only very high dimensional data might benefit from performing dimension reduction method prior to selecting new states using kNN-AS.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The effect of hyperparameters on exploration of artificial landscape. A) Artificial energy landscape with multiple local minima and one global minima. B) Decreasing the size of subsample while keeping a fixed number of neighbours (k=5) leads to more extended exploration (9 replicate). C) Altering the number of neighbours when fixing subsample size (P=0.1) has no effect on the number of visited states (9 replicates).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4. Effect of subsampling on kNN graph and state selection.** Using the whole dataset (left), the four neighbours of the datapoint (indicated with the black arrow) are from the same trajectory and the overall distribution between the neighbour vertices are skewed to 0 and 180. When selecting a subset of the data, the four neighbours for the selected datapoint are from different trajectories. This results in a uniform distribution of angles between neighbour vertices.</figcaption></figure>

The other hyperparameter, number of neighbours (k), has no clear effect on the total exploration (Figure 3C). When visualizing the states selected, we do see that a high value of k favors states at the outer boundary of the explored space, whereas a low value of k also selects states bordering “local” edges or open patches within the explored space, for example the gap
between the small and big point cloud (Figure 5). This can be compared with the UMAP algorithm
where a high number of neighbours extracts the global structure of the data while a low number
concentrates on the local structure and finer details of the data²⁸. Nonetheless, the number of
neighbors has a subtle impact on total exploration, suggesting that variations in the number of
neighbors likely only affect local exploration.
<figure><img src="image_7.png" /><figcaption>Figure 5. Differences in MD starting state selection when choosing 2, 3, 4, 5, 10, 20 or 50 neighbours. Explored states are shown in blue. Starting states for the next round of MD</figcaption></figure>
simulations are shown in orange. Last panel: Number of selected starting states not on the outer edge of the explored space.

Other hyperparameters, such as number of MD simulations per epoch, number of timesteps per simulation, and the number of epochs, are universal to any adaptive sampling strategy, and are highly dependent on the computer resources available and the system that is explored. For example, the power of adaptive sampling approach is that MD simulations can be run in parallel, and therefore, higher number of simulations per generation is better as more area can be explored within the same wall-clock time. However, the realistic number of simulations is set by the amount of resources available. Additionally, each individual trajectory must be long enough to capture new states while remaining short enough to take advantage of an adaptive sampling strategy. Lastly, the number of epochs can be determined based on a convergence criterion, for example, the extent of the area explored. For more consideration on these hyperparameters, we would refer to reviews and other adaptive sampling papers and review<sup>17,20-25,29-31</sup>. See the methods section for further details on these hyperparameters used in the experiments conducted in this paper.

## Artificial energy potentials

To evaluate the performance of kNN-AS (using either the sum of vectors or average distance metrics) in comparison to MAREAP, which is considered the state of the art in reinforcement learning-based adaptive sampling<sup>21</sup>, we first considered four simple artificial energy function test cases (Figure 6). In these test cases, we used potentials where a single particle diffuses following a Langevin dynamic. The first two potentials are a symmetric cross potential where the energy is equal across all arms and a asymmetric cross potential where there is an energy well in
one of the arms. These are key test-cases used to demonstrate the performance of MAREAP²¹. The simulations started on the horizontal arms of the energy landscape and the performance was measured by calculating the area explored by each method. The explored area is calculated as the total number of grid cells with energy level below threshold (grid size is 0.05) divided by the number of cells visited at least once. Each simulation was performed 100 times to generate error bounds. For the kNN-AS algorithm the k was set to 5 and P to 0.1.

In the symmetric cross landscape, kNN-AS (with either metric) and MAREAP explored at similar speed for the first 50 epochs (Figure 6A), but MAREAP tends to jam at the extremity of the cross, leading to only 7 out of 100 replicates exploring 90% or more of the landscape. kNN-AS on the other hand explored 90% of the landscape in 89 out of 100 replicates. Similarly with the asymmetric cross landscape, MAREAP reached 90% exploration in 1 replicate where all replicates in kNN-AS explored 90% of the landscape or more (Figure 6B). In these landscapes, there is no substantial difference observed between using the sum of vectors metric or the average distance metric.

We also investigated the versatility of the kNN-AS algorithm with more complex non-linear landscapes. The wave potential is a landscape based on a *sin* wave function where the energy is equal across the wave. The landscape is considered complex due to the quick directional changes. MAREAP explored most of the landscape without getting trapped in corners, hence the exploration curve did not plateau (Figure 6C). However, plotting state distributions reveals that simulations over-sampled the top and bottom of the wave near the starting points (Figure 6C). kNN-AS exploration was more efficient, with 60% of the landscape visited within 66 epochs versus 240 epochs on average for MAREAP. Since the energy in this system is uniform across the wave, it does not reflect real-life test cases where energy minima can be separated by high-energy
transition states. To test a more realistic energy landscape, we generated a multi-well landscape with multiple local minima and one global minimum with higher energy transitions connecting the lower energy states. Here, MAREAP often found the local minima adjacent to the starting points but failed to find the global minimum in 64 out of 100 replicates. Also, only seventeen replicates visited more than 90% of the landscape (Figure 6D). kNN-AS always found the global minimum, but still did not always explore the whole landscape as only 69 replicates visited more than 90% of the landscape. As in the simpler landscapes, the two metrics for kNN-AS show comparable performance in landscape exploration. Overall, these artificial potential functions showed that kNN-AS (regardless of the chosen metric) explores linear and complex energy landscapes efficiently and fast.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** kNN-AS efficiently explores complex artificial potential functions. A) symmetric cross potential, B) asymmetric cross potential, left plots show visualization of the energy landscape</figcaption></figure>
with stars indicating starting locations and right plot shows the area explored for kNN-AS and MAREAP over time. C) wave potential, and D) multi well potential, the top left shows a visualization of the landscape with stars indicating the starting locations. Bottom left panels show the area explored over time with kNN-AS and MAREAP. The right-hand panels show the distribution of the landscape explored by both kNN-AS (sum of vectors) and MAREAP. The color bar indicates the number of times the state was sampled.

## Adenylate Kinase

Adenylate kinase is a protein consisting of three domains and is mostly studied for its transition pathway from closed to open state. This transition can be described by the angles between the smaller domains (LID (residues 30-78) and NMP (residues 116-164)) and the larger core domain. We set up an experiment by using the open and closed crystal structures as starting points (PDB: 4AKE³² and 1AKE³³ respectively) and the angles as the features/CVs for exploration. MAREAP and kNN-AS (with either metric) were run for 400 epochs in triplicate as described in the Methods section with k set to 5 and P to 0.5 for the kNN-AS algorithm. A continuous MD simulation was also run as a baseline. MAREAP and kNN-AS were both run for an aggregated simulation time of 240 ns (6 simulations of 0.1 ns per epoch for 400 epochs), half starting from the open kinase state and half from the closed kinase conformation. To generate a similar aggregated simulation time for the continuous MD simulation, two simulations are started from each starting point for 120 ns.

After 240 ns, kNN-AS had explored on average 48% and 52% of the space (sum of vectors and average distance respectively), MAREAP 36% and full MD 20% (Figure 7A). This shows that
both MAREAP and kNN-AS explored significantly more area than the baseline MD simulation. Interestingly, both MAREAP and kNN-AS (with either metric) retrieved an alternate crystal structure (PDB: 2AK3³⁴) that was missed by the regular MD simulation (Figure 7B). Additionally, kNN-AS using the average distance metric results in a more densely filled space, while the sum of vectors leads to a more sparse exploration (Figure 7B). To evaluate whether kNN-AS pushes the simulations to unrealistic conformations, we visualized the structures at the edges of the NMP-LID angle subspace (Figure S2) and did not observe signs of irregularities such as alpha distorted helices/beta sheets or unfolding of protein regions.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Results on Adenylate Kinase conformational landscape. A) Percentage of area explored vs epochs. B) Visualization of area explored in LID/NMP angle space for each method, when starting from open (PDB: 4AKE³²).) and closed (PDB: 1AKE³³) conformations. Experimentally observed alternate conformations are shown on the landscape (PDB: 2AK3³⁴ and 1DVR³⁵).</figcaption></figure>

Since simulations have been restarted from the edges, the distribution from the raw data contains bias and noise. To assess the true distributions and kinetics from data generated using
adaptive sampling methods, MSM models can be constructed from the off-equilibrium trajectories,
as no biasing force was applied during the simulations¹⁸. These models reveal the stationery and
transition probabilities of states and can reweight trajectories to correct for the bias introduced by
over-sampling less populated regions. We ran kNN-AS (sum of vectors) with longer simulation
time per trajectory and clustered the frames based on distances between the three domains (LID-
NMP-CORE). The MSM model showed convergence for the slowest process for higher lagtimes
(Figure S3A). After reweighting the trajectories, we plotted the free energy landscape on the LID
and NMP angles. We observe that apo/ligand-free AdK is mostly in the open conformation (Figure
S3B) which is similar to what is observed experimentally²³. Also, two transition states are observed
with either one of the domains closed and the other one open, where the state with closed-LID and
open-NMP sample more frequently which has been demonstrated previously³⁶,³⁷. The slowest
process corresponds to the closing of the NMP domain (Figure S3C) which has been described as
the rate-limiting step for conformational transitions before³⁸.

Discussion

MD simulations of large biomolecular systems, such as multi-protein complexes, often
require impractical amounts of resources to reach equilibrium. Therefore, optimized strategies to
sample the free-energy landscape aim to accelerate this process by either running parallel
simulations with carefully chosen starting states (adaptive sampling)²⁰,²¹,³⁰,³¹, or by introducing a
bias to push the simulations towards underexplored areas or other regions of interest²⁰,²¹,³⁰,³¹. The
latter requires reweighing of the energy landscape¹⁴ while adaptive sampling methods are run
under normal conditions and can be easily analyzed using Markov State Models for underlying
dynamics¹⁹. Using artificial neural networks to learn quantum-mechanical force fields or rapidly
generate ab initio protein equilibrium ensembles are also novel and promising avenues to
accurately and efficiently sample conformational states³⁹⁻⁴². In this application note, we present a
very simple adaptive sampling strategy where a kNN-graph detects states on the edge of the
conformational landscape and determines which areas have been explored.

Since kNN-AS is based on the kNN graph, it may suffer from the curse of dimensionality
when using high-dimensional CV space⁴³. Also, kNN-AS efficiently explores the space defined by
a set of CVs, but it does not learn which CVs are most important for sampling, unlike MAREAP.
Therefore, a future direction could be to pair kNN-AS with dimension reduction techniques, such
as principal component analysis⁴⁴, variational mixture of posteriors⁴⁵ or variational
autoencoders⁴⁶,⁴⁷ to extract a lower CVs space for more efficient exploration.

An important hyperparameter is the size of subsample used in constructing the kNN-graph.
Using a subset introduces a level of randomness in the process allowing for heavily explored areas
to be explored again if only a small set is sampled from the region. This can be advantageous for
landscapes with high energy transition states, as MD simulations need to be restarted multiple
times at the same edge before a transition is observed. Importantly, the size of subsample is
dependent on the number of frames saved during the simulations. Shorter intervals cause frames to
be too close in space, resulting in near neighbours always coming from the same trajectory and
limiting efficient exploration across trajectories. Therefore, it is recommended to either use small
data samples when generating k-NN graphs, or keeping long intervals between saved frames, or
both. For the test cases shown here, subsamples of 10% and 50% of the data were used for the
artificial energy potentials and the adenylate kinase cases respectively. The lower number of
frames required for kNN-AS to perform well has benefits for simulating large biomolecular system
as less data needs to be stored in memory.
The number of neighbours in kNN is a hyperparameter that influences state selection on local versus global edges (Figure 4), but has less of an effect on the total exploration. We decided to use five nearest neighbours, since this value of k showed good performance on the artificial landscapes, but future research could further interrogate how focus on local versus global exploration differently affects results of MD simulations. Since kNN-AS only requires hyperparameter tuning for 2 parameters, k and P (aside from regular adaptive sampling hyperparameters), a few combinations can be tested on shorter runs for a given system before launching longer simulations.

In this paper, we explored two metrics for kNN-AS to identify boundary points, namely the sum of vectors and average distance. Both the sum of vectors and average distance methods identify edges in the landscape. However, the former also highlights corners, which may represent potential transition pathways. In contrast, the average distance approach emphasizes sparsely sampled regions. Only in the AdK test case, a marginal improvement is observed for the average distance metric, which could be due to filling in more of the sparse areas. An interesting future research idea would be to explore switching between the two metrics to optimally explore path transitions as well as low sampled areas.

Overall, we have shown that kNN-AS is a simple and short algorithm that efficiently guides the exploration of complex conformational landscapes with minimal hyperparameter tuning. This simple code can easily be inserted within more detailed algorithms, can be augmented with complementary techniques such as scores re-weighting for directed searches, or combined with deep learning approaches such as VAEs.
# Methods

## Langevin Dynamics

Four artificial energy potential functions were chosen to compare the performance of kNN-AS with MAREAP. The first two potentials (symmetric cross and asymmetric cross) have been previously described and tested for the MAREAP algorithm²¹. The other two energy potentials were developed to mimic more complex energy landscapes (wave and multi well potentials). The energy potentials are defined as:

### Symmetric cross potential²¹

$$V(x, y) = -25 * e^{-0.5 * (|x - 0.2|^2 / 0.008 + |y - 1|^2 / 0.008)} - 25 * e^{-0.5 * (|x - 1|^2 / 0.008 + |y - 0.2|^2 / 0.008)} - 25 * e^{-0.5 * (|x - 1.8|^2 / 0.008 + |y - 1|^2 / 0.008)} - 25 * e^{-0.5 * (|x - 2|^2 / 0.008 + |y - 2|^2 / 0.008)}$$

### Asymmetric cross potential²¹

$$V(x, y) = -80 * e^{-0.5 * (|x - 0.2|^2 / 0.008 + |y - 1|^2 / 0.008)} - 25 * e^{-0.5 * (|x - 1|^2 / 0.008 + |y - 0.2|^2 / 0.008)} - 25 * e^{-0.5 * (|x - 1.8|^2 / 0.008 + |y - 1|^2 / 0.008)} - 25 * e^{-0.5 * (|x - 2|^2 / 0.008 + |y - 2|^2 / 0.008)}$$

### Wave potential

$$V(x, y) = -80 * e^{-[20 * (1 - \tanh(8 - |x*3|)) + 0.5 * (|y*2 - \sin(2 * \pi * |x*3|) / 4| / 0.8)^2]}$$

### Multiwell potential

$$V(x, y) = -60 * e^{-1 * (|x - (-1.1)|^2 - 1.5 * |y - 0|)^2} \pm 60 * e^{-2 * (|x - 0.8|^2 - 2 * |y - (-0.3)|)^2} \pm 60 * e^{-6 * (|x - (-1.1)|^2 - 3.5 * |y - (-1.2)|)^2} \pm 60 * e^{-5 * (|x - (1.55)|^2 - 5)^2}$$

All simulations were run with OpenMM 7.7⁴⁸ and a Langevin integrator. The timestep was set to 2 fs, the friction coefficient to 1 ps⁻¹ and the temperature to 300K. The mass of the single
particle was 100 Da. For both protocols, MAREAP and kNN-AS, 400 epochs were run with 6 MD
simulations for each epoch (each running for 500 steps) and a frame was saved every 25 steps. Only
the symmetric cross potential was run for 200 epochs. For each test case, 100 replicates were
performed to determine error bounds. The starting locations are: symmetric cross [[0.8, 1], [1.2,
1]], asymmetric cross [[0.2, 1], [1.8, 1]], wave [[-1.25, 0.5], [1.25, -0.5]], and multi-well [[1.6, 0.4],
[-1.2, -1.2]].

For the hyperparameter sensitivity study, the multi well landscape was chosen, and the
kNN-AS (sum of vectors) protocol was run 400 epochs with ranging k [2,4,5,10] and ranging P
[0.1, 0.2, 0.6, 1.0]. P was set to 0.1 for the ranging k experiments and k was set to 5 for the ranging
P studies.

Hyperparameters for MAREAP were similar to the parameters previously described²¹ with
10⁴ frames used for κ-means clustering (with 5 restarts for random cluster initialization). The
number of clusters was determined using the CLUST⁴⁹ calculation by using d=2, γ=0.6 and b=3⁻⁴.
The maximum weight change was set to δ=0.02, and collective rewards and fraction stakes were
used. For kNN-AS, 10% of the data was subsampled for kNN graph generation and number of
neighbours was set to 5.

All-atom dynamics

MD simulations for Adenylate Kinase were also performed with OpenMM 7.7⁴⁸ with a
Langevin integrator, temperature of 300K, time step of 2 fs and friction coefficient of 1 ps⁻¹. All
bonds with hydrogens were constrained and the pressure was held at 1 bar using the Monte Carlo
barostat. The starting states are the open and closed conformations (PDB: 4AKE<sup>32</sup> and 1AKE<sup>33</sup> respectively) and were prepared using the Amber 14 forcefield<sup>50</sup> and TIP3-pfb<sup>51</sup> water model. 40000 water molecules were added to the systems to obtain a similar number of atoms in both systems. The ionic strength was set to 0.1M. Then the system was energy-minimized and equilibrated for 100 steps. The LID angle was calculated for the vectors from residues 115-125 to 125-153 and 179-185. The NMP angle was calculated for the vectors from residues 90-99 to 115-125 and 35-55.

MAREAP and kNN-AS were both run for 400 epochs with 6 MD simulations per epoch. Each MD simulation consisted of 50,000 steps and positions were saved every 3500 steps. To generate error-bounds, each method was run in triplicate. MAREAP was initialized using similar hyperparameters as the artificial energy function along with using 12 least-count clusters. kNN-AS was run with 5 neighbours for kNN graph generation and 50% of the data was used. Regular MD simulations were performed starting from the two starting points each for 120ns to get a total of 240ns aggregated simulation time, similar to MAREAP and kNN-AS.

## MSM model generation

The experiment for AdK was repeated with longer simulation time (20ns per simulation) for 6 simulations per generation for a total of 40 generations. A frame was saved every 10ps. We clustered the trajectories based on the Cα-distances between LID-CORE, NMP-CORE, and LID-NMP domains in 35 clusters using the K-means++ algorithm from deeptime using 20 iterations. The final MSM model was generated using a lagtime of 4 ns (400 steps).
## Supporting information

Supplementary figures showing the effect of increasing dataset or dimension size on computation time, configurations of Adenylate Kinase at the edge of the LID-NMP angle subspace, and a MSM model estimated on the data generated by kNN-AS.

## Data availability

The scripts to reproduce the results presented and the script containing the kNN-AS algorithm is available at: https://github.com/ERovers/kNN-AS. Raw results are available on request.

## Acknowledgements

This research was enabled in part by support provided by Compute Ontario (https://www.computeontario.ca/) and the Digital Research Alliance of Canada (alliancecan.ca). M.S. and E.R. gratefully acknowledges support from NSERC [Grant RGPIN-2019-04416], CQDM (Quantum Leap-176), and MITACS accelerate (IT13051). The Structural Genomics Consortium is a registered charity (no: 1097737) that receives funds from Bayer AG, Boehringer Ingelheim, Bristol Myers Squibb, Genentech, Genome Canada through Ontario Genomics Institute [OGI-196], EU/EFPIA/OICR/McGill/KTH/Diamond Innovative Medicines Initiative 2 Joint Undertaking [EUbOPEN grant 875510], Janssen, Merck KGaA (aka EMD in Canada and US), Pfizer, and Takeda. C.M. is supported in part by the Province of Ontario, the Government of Canada through CIFAR, and companies sponsoring the Vector Institute and the support of the
Natural Sciences and Engineering Research Council of Canada (NSERC), RGPIN-2021-03445.

A.T. is supported by a Vanier Fellowship from NSERC.
References

(1) Bhayye, S. S.; Roy, K.; Saha, A. Molecular Dynamics Simulation Study Reveals Polar Nature of Pathogenic Mutations Responsible for Stabilizing Active Conformation of Kinase Domain in Leucine-Rich Repeat Kinase II. *Struct Chem* **2018**, *29* (3), 657–666.
https://doi.org/10.1007/s11224-017-1059-z

(2) Khare, S. D.; Ding, F.; Dokholyan, N. V. Folding of Cu, Zn Superoxide Dismutase and Familial Amyotrophic Lateral Sclerosis. arXiv June 30, 2003.
https://doi.org/10.48550/arXiv.physics/0306200

(3) Khazamipour, A.; Gholampour-Faroji, N.; Zeraati, T.; Vakilian, F.; Haddad-Mashadrizeh, A.; Ghayour Mobarhan, M.; Pasdar, A. A Novel Causative Functional Mutation in GATA6 Gene Is Responsible for Familial Dilated Cardiomyopathy as Supported by in Silico Functional Analysis. *Sci Rep* **2022**, *12* (1), 13752. https://doi.org/10.1038/s41598-022-13993-6

(4) Odera, M.; Furuta, T.; Sohma, Y.; Sakurai, M. Molecular Dynamics Simulation Study on the Structural Instability of the Most Common Cystic Fibrosis-Associated Mutant ΔF508-CFTR. *Biophysics and Physicobiology* **2018**, *15*, 33–44.
https://doi.org/10.2142/biophysico.15.0_33

(5) Wang, F.; Orioli, S.; Ianeselli, A.; Spagnolli, G.; Beccara, S. a; Gershenson, A.; Faccioli, P.; Wintrode, P. L. All-Atom Simulations Reveal How Single Point Mutations Promote Serpin Misfolding. *Biophysical Journal* **2018**, *114* (9), 2083–2094.
https://doi.org/10.1016/j.bpj.2018.03.027

(6) Wilson, C. J.; Chang, M.; Karttunen, M.; Choy, W.-Y. KEAP1 Cancer Mutants: A Large-Scale Molecular Dynamics Study of Protein Stability. *International Journal of Molecular Sciences* **2021**, *22* (10), 5408. https://doi.org/10.3390/ijms22105408

(7) Comitani, F.; Gervasio, F. L. Exploring Cryptic Pockets Formation in Targets of Pharmaceutical Interest with SWISH. *J. Chem. Theory Comput.* **2018**, *14* (6), 3321–3331.
https://doi.org/10.1021/acs.jctc.8b00263

(8) Kuzmanic, A.; Bowman, G. R.; Juarez-Jimenez, J.; Michel, J.; Gervasio, F. L. Investigating Cryptic Binding Sites by Molecular Dynamics Simulations. *Acc. Chem. Res.* **2020**, *53* (3), 654–661. https://doi.org/10.1021/acs.accounts.9b00613

(9) Meller, A.; Ward, M.; Borowsky, J.; Kshirsagar, M.; Lotthammer, J. M.; Oviedo, F.; Ferres, J. L.; Bowman, G. R. Predicting Locations of Cryptic Pockets from Single Protein Structures Using the PocketMiner Graph Neural Network. *Nat Commun* **2023**, *14* (1), 1177. https://doi.org/10.1038/s41467-023-36699-3

(10) Raich, L.; Meier, K.; Günther, J.; Christ, C. D.; Noé, F.; Olsson, S. Discovery of a Hidden Transient State in All Bromodomain Families. *Proceedings of the National Academy of Sciences* **2021**, *118* (4), e2017427118. https://doi.org/10.1073/pnas.2017427118

(11) Sztain, T.; Amaro, R.; McCammon, J. A. Elucidation of Cryptic and Allosteric Pockets within the SARS-CoV-2 Main Protease. *J. Chem. Inf. Model.* **2021**, *61* (7), 3495–3501.
https://doi.org/10.1021/acs.jcim.1c00140

(12) Bowman, G. R.; Geissler, P. L. Equilibrium Fluctuations of a Single Folded Protein Reveal a Multitude of Potential Cryptic Allosteric Sites. *Proceedings of the National Academy of Sciences* **2012**, *109* (29), 11681–11686. https://doi.org/10.1073/pnas.1209309109
(13) Knoverek, C. R.; Amarasinghe, G. K.; Bowman, G. R. Advanced Methods for Accessing Protein Shape-Shifting Present New Therapeutic Opportunities. *Trends in Biochemical Sciences* **2019**, *44* (4), 351–364. https://doi.org/10.1016/j.tibs.2018.11.007.

(14) Hénin, J.; Lelièvre, T.; Shirts, M. R.; Valsson, O.; Delemotte, L. Enhanced Sampling Methods for Molecular Dynamics Simulations. *LiveCoMS* **2022**, *4* (1). https://doi.org/10.33011/livecoms.4.1.1583.

(15) Ray, D.; Ansari, N.; Rizzi, V.; Invernizzi, M.; Parrinello, M. Rare Event Kinetics from Adaptive Bias Enhanced Sampling. *J. Chem. Theory Comput.* **2022**, *18* (11), 6500–6509. https://doi.org/10.1021/acs.jctc.2c00806.

(16) Lan, H.; Hsia, O.; Nakasone, M.; Wijaya, A.; Ciulli, A.; Butts, C.; Aggarwal, V. Computation, Synthesis and NMR Investigations of PROTAC Linker Conformation, Chamaleonicity and Their Impacts on the Mode of Action. *ChemRxiv* October 23, 2024. https://doi.org/10.26434/chemrxiv-2024-5dbnb.

(17) Kleiman, D. E.; Nadeem, H.; Shukla, D. Adaptive Sampling Methods for Molecular Dynamics in the Era of Machine Learning. *J. Phys. Chem. B* **2023**, *127* (50), 10669–10681. https://doi.org/10.1021/acs.jpcb.3c04843.

(18) Noé, F.; Schütte, C.; Vanden-Eijnden, E.; Reich, L.; Weikl, T. R. Constructing the Equilibrium Ensemble of Folding Pathways from Short Off-Equilibrium Simulations. *Proceedings of the National Academy of Sciences* **2009**, *106* (45), 19011–19016. https://doi.org/10.1073/pnas.0905466106.

(19) *An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation*; Bowman, G. R., Pande, V. S., Noé, F., Eds.; Advances in Experimental Medicine and Biology; Springer Netherlands: Dordrecht, 2014; Vol. 797. https://doi.org/10.1007/978-94-007-7606-7.

(20) Pérez, A.; Herrera-Nieto, P.; Doerr, S.; De Fabritiis, G. Adaptive Bandit: A Multi-Armed Bandit Framework for Adaptive Sampling in Molecular Simulations. *J. Chem. Theory Comput.* **2020**, *16* (7), 4685–4693. https://doi.org/10.1021/acs.jctc.0c00205.

(21) Kleiman, D. E.; Shukla, D. Multiagent Reinforcement Learning-Based Adaptive Sampling for Conformational Dynamics of Proteins. *J Chem Theory Comput* **2022**, *18* (9), 5422–5434. https://doi.org/10.1021/acs.jctc.2c00683.

(22) Inizan, T. J.; Célerse, F.; Adjoua, O.; Ahdab, D. E.; Jolly, L.-H.; Liu, C.; Ren, P.; Montes, M.; Lagarde, N.; Lagardère, L.; Monmarché, P.; Piquemal, J.-P. High-Resolution Mining of the SARS-CoV-2 Main Protease Conformational Space: Supercomputer-Driven Unsupervised Adaptive Sampling. *Chem. Sci.* **2021**, *12* (13), 4889–4907. https://doi.org/10.1039/D1SC00145K.

(23) Chiavazzo, E.; Covino, R.; Coifman, R. R.; Gear, C. W.; Georgiou, A. S.; Hummer, G.; Kevrekidis, I. G. Intrinsic Map Dynamics Exploration for Uncharted Effective Free-Energy Landscapes. *Proceedings of the National Academy of Sciences* **2017**, *114* (28), E5494–E5503. https://doi.org/10.1073/pnas.1621481114.

(24) Harada, R.; Nakamura, T.; Takano, Y.; Shigeta, Y. Protein Folding Pathways Extracted by OFLOOD: Outlier FLOODing Method. *Journal of Computational Chemistry* **2015**, *36* (2), 97–102. https://doi.org/10.1002/jcc.23773.

(25) Harada, R.; Shigeta, Y. Hybrid Cascade-Type Molecular Dynamics with a Markov State Model for Efficient Free Energy Calculations. *J. Chem. Theory Comput.* **2019**, *15* (1), 680–687. https://doi.org/10.1021/acs.jctc.8b00802.
(26) Le, M. Q.; Taylor, D. Persistent Homology with K-Nearest-Neighbor Filtrations Reveals Topological Convergence of PageRank. arXiv June 9, 2022.
https://doi.org/10.48550/arXiv.2206.04725.

(27) Carli, M.; Rodriguez, A.; Laio, A.; Glielmo, A. Density Estimation via Binless Multidimensional Integration. arXiv July 14, 2024.
https://doi.org/10.48550/arXiv.2407.08094.

(28) McInnes, L.; Healy, J.; Melville, J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv September 18, 2020.
https://doi.org/10.48550/arXiv.1802.03426.

(29) Betz, R. M.; Dror, R. O. How Effectively Can Adaptive Sampling Methods Capture Spontaneous Ligand Binding? *J Chem Theory Comput* **2019**, *15* (3), 2053–2063.
https://doi.org/10.1021/acs.jctc.8b00913.

(30) Zimmerman, M. I.; Bowman, G. R. FAST Conformational Searches by Balancing Exploration/Exploitation Trade-Offs. *J Chem Theory Comput* **2015**, *11* (12), 5747–5757.
https://doi.org/10.1021/acs.jctc.5b00737.

(31) Shamsi, Z.; Cheng, K. J.; Shukla, D. Reinforcement Learning Based Adaptive Sampling: REAPing Rewards by Exploring Protein Conformational Landscapes. *J Phys Chem B* **2018**, *122* (35), 8386–8395. https://doi.org/10.1021/acs.jpcb.8b06521.

(32) Müller, C.; Schlauderer, G.; Reinstein, J.; Schulz, G. Adenylate Kinase Motions during Catalysis: An Energetic Counterweight Balancing Substrate Binding. *Structure* **1996**, *4* (2), 147–156. https://doi.org/10.1016/S0969-2126(96)00018-4.

(33) Müller, C. W.; Schulz, G. E. Structure of the Complex between Adenylate Kinase from Escherichia Coli and the Inhibitor Ap5A Refined at 1.9 Å Resolution. A Model for a Catalytic Transition State. *J Mol Biol* **1992**, *224* (1), 159–177. https://doi.org/10.1016/0022-2836(92)90582-5.

(34) Diederichs, K.; Schulz, G. E. The Refined Structure of the Complex between Adenylate Kinase from Beef Heart Mitochondrial Matrix and Its Substrate AMP at 1.85 Å Resolution. *Journal of Molecular Biology* **1991**, *217* (3), 541–549. https://doi.org/10.1016/0022-2836(91)90756-V.

(35) Schlauderer, G. J.; Proba, K.; Schulz, G. E. Structure of a Mutant Adenylate Kinase Ligated with an ATP-Analogue Showing Domain Closure Over ATP. *Journal of Molecular Biology* **1996**, *256* (2), 223–227. https://doi.org/10.1006/jmbi.1996.0080.

(36) Lee, J.; Joo, K.; Brooks, B. R.; Lee, J. The Atomistic Mechanism of Conformational Transition of Adenylate Kinase Investigated by Lorentzian Structure-Based Potential. *J. Chem. Theory Comput.* **2015**, *11* (7), 3211–3224. https://doi.org/10.1021/acs.jctc.5b00268.

(37) Zheng, Y.; Cui, Q. Multiple Pathways and Time Scales for Conformational Transitions in Apo-Adenylate Kinase. *J. Chem. Theory Comput.* **2018**, *14* (3), 1716–1726.
https://doi.org/10.1021/acs.jctc.7b01064.

(38) Whitford, P. C.; Miyashita, O.; Levy, Y.; Onuchic, J. N. Conformational Transitions of Adenylate Kinase: Switching by Cracking. *J Mol Biol* **2007**, *366* (5), 1661–1671.
https://doi.org/10.1016/j.jmb.2006.11.085.

(39) Diez, J. V.; Schreiner, M.; Engkvist, O.; Olsson, S. Boltzmann Priors for Implicit Transfer Operators. arXiv October 14, 2024. https://doi.org/10.48550/arXiv.2410.10605.

(40) Jing, B.; Stärk, H.; Jaakkola, T.; Berger, B. Generative Modeling of Molecular Dynamics Trajectories. arXiv September 26, 2024. https://doi.org/10.48550/arXiv.2409.17808.
(41) Noé, F.; Olsson, S.; Köhler, J.; Wu, H. Boltzmann Generators: Sampling Equilibrium States of Many-Body Systems with Deep Learning. *Science* **2019**, *365* (6457), eaaw1147.
https://doi.org/10.1126/science.aaw1147.

(42) Unke, O. T.; Stöhr, M.; Ganscha, S.; Unterthiner, T.; Maennel, H.; Kashubin, S.; Ahlin, D.; Gastegger, M.; Medrano Sandonas, L.; Berryman, J. T.; Tkatchenko, A.; Müller, K.-R. Biomolecular Dynamics with Machine-Learned Quantum-Mechanical Force Fields Trained on Diverse Chemical Fragments. *Science Advances* **2024**, *10* (14), eadn4397.
https://doi.org/10.1126/sciadv.adn4397.

(43) Bellman, R. INTRODUCTION. In *Adaptive Control Processes: A Guided Tour*; Princeton University Press, 1961; pp 3–12.

(44) Pearson, K. LIII. On Lines and Planes of Closest Fit to Systems of Points in Space. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science* **1901**, *2* (11), 559–572. https://doi.org/10.1080/14786440109462720.

(45) Tomczak, J. M.; Welling, M. VAE with a VampPrior. arXiv February 26, 2018. https://doi.org/10.48550/arXiv.1705.07120.

(46) Kingma, D. P.; Welling, M. Auto-Encoding Variational Bayes. arXiv December 20, 2013. https://doi.org/10.48550/arXiv.1312.6114.

(47) Kingma, D. P.; Welling, M. Auto-Encoding Variational Bayes. arXiv December 10, 2022. https://doi.org/10.48550/arXiv.1312.6114.

(48) Eastman, P.; Swails, J.; Chodera, J. D.; McGibbon, R. T.; Zhao, Y.; Beauchamp, K. A.; Wang, L. P.; Simmonett, A. C.; Harrigan, M. P.; Stern, C. D.; Wiewiora, R. P.; Brooks, B. R.; Pande, V. S. OpenMM 7: Rapid Development of High Performance Algorithms for Molecular Dynamics. *PLoS computational biology* **2017**, *13* (7).
https://doi.org/10.1371/JOURNAL.PCBI.1005659.

(49) Buenfil, J.; Koelle, S. J.; Meila, M. Tangent Space Least Adaptive Clustering; 2021.

(50) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. *J. Chem. Theory Comput.* **2015**, *11* (8), 3696–3713.
https://doi.org/10.1021/acs.jctc.5b00255.

(51) Wang, L.-P.; Martinez, T. J.; Pande, V. S. Building Force Fields: An Automatic, Systematic, and Reproducible Approach. *J. Phys. Chem. Lett.* **2014**, *5* (11), 1885–1891.
https://doi.org/10.1021/jz500737m.

For Table of Contents Only

<figure><img src="image_10.png" /></figure>