# Automating Staged Rollout with Reinforcement Learning

Shadow Pritchard  University of TulsaTulsaOklahomaUSA<EMAIL_ADDRESS>,
Vidhyashree Nagaraju  University of TulsaTulsaOklahomaUSA vidhyashree-
<EMAIL_ADDRESS>and  Lance Fiondella  University of Massachusetts
DartmouthDartmouthMassachusettsUSA<EMAIL_ADDRESS>
(2022)

###### Abstract.

Staged rollout is a strategy of incrementally releasing software updates to
portions of the user population in order to accelerate defect discovery
without incurring catastrophic outcomes such as system wide outages. Some past
studies have examined how to quantify and automate staged rollout, but stop
short of simultaneously considering multiple product or process metrics
explicitly. This paper demonstrates the potential to automate staged rollout
with multi-objective reinforcement learning in order to dynamically balance
stakeholder needs such as time to deliver new features and downtime incurred
by failures due to latent defects.

DevOps, Staged Rollout, Reinforcement Learning, Software Reliability

††journalyear: 2022††copyright: acmcopyright††conference: New Ideas and
Emerging Results ; May 21–29, 2022; Pittsburgh, PA, USA††booktitle: New Ideas
and Emerging Results (ICSE-NIER’22), May 21–29, 2022, Pittsburgh, PA,
USA††price: 15.00††doi: 10.1145/3510455.3512782††isbn:
978-1-4503-9224-2/22/05††ccs: Software and its engineering Software testing
and debugging

##  1\. Introduction

Researchers recognize the importance of software performance engineering and
application performance management (Brunnert et al., 2015) during the
development and operation phases, including interoperability between different
tools and techniques. However, it is also important to consider process and
product quality metrics in an integrated manner as well as any tradeoffs
imposed. For example, staged rollout (Zhao et al., 2018) seeks to accelerate
deployment of new features while preserving safety, yet this strategy of
introducing an update to a subset of users may pose tradeoffs such as the time
to deliver a new feature and the downtime experienced by users. Techniques to
automate the staged rollout process could reduce the level of expertise
required to balance these stakeholder needs.

Tarvo et al. (Tarvo et al., 2015) defined and collected metrics on software
updates, but required the user to make release decisions. Zhao et al. (Zhao et
al., 2018) introduced a framework to release software to subsets of users
based on time, power, and risk-based scheduling, which rely on predetermined
thresholds or time windows. The thesis of Velayutham (Velayutham, 2021)
treated staged rollout as a time series problem with autoregressive integrated
moving average and long short-term memory to automate decision making, but
tradeoffs were not explicitly considered.

Potential reinforcement learning approaches to model staged rollout include
non-stationary (Khetarpal et al., 2020) methods, which encompass time-varying
rewards. Abdallah et al. (Abdallah and Kaisers, 2016) introduced repeated
Q-learning to update reward tables more frequently for seldom visited states
in order to improve performance under non-stationary conditions. Multi-
objective reinforcement learning (Liu et al., 2014) attempts to balance
rewards associated with multiple objectives through a weighted, linear reward
(Krass, 1990). Additional methods to balance multiple objectives include
geometric steering (Vamplew et al., 2017) or hyper volume action selection
(Van Moffaert et al., 2013).

This paper presents a model of the staged rollout problem as a non-stationary
Markov decision process (MDP). Multi-objective Q-learning with upper
confidence bound exploration is applied to solve the MDP. The approach is
demonstrated on a data set from the software reliability engineering
literature, illustrating reinforcement learning as a promising approach to
make dynamic decisions during staged rollout.

The remainder of the paper is organized as follows: Section 2 formulates a
state model of the staged rollout process, discussing tradeoffs and reward
function modeling. Section 3 summarizes Q-learning with upper confidence bound
as well as naive policy enumeration. Section 4 develops metrics to assess
policies identified by reinforcement learning. Section 5 presents preliminary
results. Section 6 describes plans for future research.

##  2\. Staged Rollout Modeling and Tradeoff Assessment

This section presents a state model of the staged rollout problem that can be
interpreted as a Markov decision process (MDP), where a multi-objective
reinforcement learning agent makes decisions to balance two primary factors,
including (i) downtime and (ii) delivery time according to stakeholder
preference. For example, safety critical software industries place greater
emphasis on minimizing failures, whereas many other industries prioritize
minimizing delivery time to get a product to market quickly.

###  2.1. Staged Rollout Model

Staged rollout of software (Zhao et al., 2018) has been recognized as a
strategy to field new functionality on an ongoing basis without incurring
failures that induce system outages, widespread unavailability of services,
economic losses, and user dissatisfaction. The rationale for staged rollout is
to publish an updated software possessing new functionality for use by a
subset of the user base to avoid the major problems described above, but also
to accelerate the discovery of defects. The development team then attempts to
correct the source of the problem and begins the process of staged rollout
anew.

Figure 1 shows a simple state diagram of a staged rollout process.

Figure 1. State diagram of staged rollout

State $Dev$ represents the development state, where software is tested by an
internal team. An elementary model of the traditional approach to software
updates simply transitions to the $Ops$ state, once the software is deemed
satisfactory with respect to functional requirements, reliability, and other
desired attributes, where the $Ops$ state exposes the software to the entire
base of $n_{Ops}$ users. Staged rollout, instead, transitions from the $Dev$
state to state $i_{1}$, which represents the first stage of staged rollout,
where the software is published for use by $p_{i_{1}}$ percent or
$n_{i_{1}}=p_{i_{1}}\times n_{Ops}$ of the user base. Multiple stages of
staged rollout between development and full deployment are also possible. This
more general case possesses $m$ intermediate staged rollout states, where it
may be reasonable to assume that each transition from $i_{j}$ to $i_{j+1}$
increases the fraction of the user base exposed to the software such that
$p_{i_{j}}<p_{i_{j+1}}$ and $n_{i_{j}}<n_{i_{j+1}}$. Failure in any state
transitions to the $Dev$ state, where root cause analysis and defect removal
are attempted.

The state model described in Figure 1 enables explicit consideration of the
tradeoffs between downtime and delivery time. Downtime is determined by the
state in which the failure occurs and is proportional to the fraction of the
user base $n_{i_{j}}$ and mean time to resolve ($MTTR$). Thus, failure in the
$j^{th}$ state of staged rollout ($i_{j}$) contributes less to downtime than
failure in the state $i_{j+1}$. However, the defect exposure rate in the
$j^{th}$ state of staged rollout is also less than the defect exposure rate in
the $(j+1)^{th}$ state, meaning that the downtime experienced, and the time
required to discover and remove all defects are modeled as competing
constraints. Thus, minimizing downtime by remaining in the $Dev$ state until
all defects have been detected and removed will likely delay delivery time.
Similarly, unrestrained transition to the Ops state immediately after each
defect is resolved is likely to exacerbate downtime. Therefore, it may not be
possible to simultaneously minimize downtime and delivery time, posing a
multi-objective problem. Moreover, organizations implementing staged rollout
may express different levels of tolerance for these undesirable outcomes.
Subsequently, it is unlikely that a single optimal policy or ”one size fits
all” approach to staged rollout exists. Instead, it is necessary to select
transition times $t_{j,j+1}=\frac{1}{\lambda_{j,j+1}}$ that balance downtime
and delivery time in a manner that is satisfactory to the customer.
Intuitively, a high defect discovery rate in the $Dev$ state is likely to
indicate that additional defects will be uncovered. Hence, staged rollout
should not be performed because it would risk greater downtime. Therefore, the
problem is to select numerical values of transition rates
$\lambda_{j,j+1}=\frac{1}{t_{j,j+1}}$ that achieve the desired balance between
downtime and delivery time. In the case where these rates vary as a function
of time, the problem is said to be non-stationary.

Since the transitions $\lambda_{Dev,i_{1}}$ and $\lambda_{i_{1},Ops}$ of
Figure 1 are decisions, the staged rollout state diagram is a Markov decision
process to which reinforcement learning can be applied in order to solve for a
range of policies to take actions (make transition decisions) from a given
state. Traditional software reliability models (Farr, 1996) assume that the
rate of defect discovery is nonhomogeneous. For example, if the rate of defect
discovery is high in the early stages of testing, then the optimal action will
be to stay in the $Dev$ state, where the penalty for failure is lowest.
However, later in testing when fewer undiscovered defects remain, the optimal
action will be to transition to state $i_{1}$ or $Ops$ in order to reduce
delivery time. Therefore, the model is non-stationary because the rate at
which defects are discovered changes as testing progresses. Thus, this change
in the defect detection rate will also cause the reward to change and the
optimal policy for a given state will change over time.

###  2.2. Process Tradeoffs and Reward Modeling

This section defines quantitative metrics for delivery time and downtime in
terms of the state model shown in Figure 1 as well as how these metrics are
incorporated into a multi-objective reward model suitable for application of
reinforcement learning.

####  2.2.1. Delivery Time

To model the impact of staged rollout on delivery time, we assume that one
unit of time in the $Dev$ state advances the timeline by one unit, whereas
time in the staged rollout and $Ops$ state accelerate the rate at which time
advances proportional to the percentage of the user base. Therefore, staged
rollout may be regarded as a modern form of accelerated life testing (Nelson,
2009) for software. For example, if the complete user base is composed of
$n_{Ops}=10,000$ users and staged rollout exposes new functionality to
$p_{i_{1}}=0.1$ or $10\%$ of the user base, then $n_{i_{1}}=1000$. Similarly,
if $n_{Dev}=50$, then a simple method to compute the acceleration factor in
each state of staged rollout is the ratio between the number of users in a
state over the baseline in the $Dev$ state such that the acceleration factor
in the staged rollout and $Ops$ states are $\frac{n_{i_{1}}}{n_{Dev}}=20$ and
$\frac{n_{Ops}}{n_{Dev}}=200$ respectively. This simplifying assumption can be
improved, since testers are familiar with the functionality and intentionally
stress the program to expose defects. Modeling these ratios is a research
question that requires staged rollout data. Nevertheless, the simplifying
assumptions made here enable a quantitative framework upon which to improve.

The preliminary assumption of linear acceleration factors described above
provides a concrete starting point to measure the cumulative time required to
reach a release decision. Specifically, delivery time may be defined as the
time to discover and eliminate defects associated with a test set plus the
time to transition from the $Dev$ to $Ops$ state
($t_{Dev,i_{1}}+t_{t_{1},Ops}$) because no additional defects are discovered
under the simplifying assumption that the final defect is resolved
immediately. Modeling advances that explicitly consider the time between
defect discovery and resolution (Lo and Huang, 2006; Nafreen et al., 2020) can
further enhance the realism of the staged rollout deployment model.

####  2.2.2. Downtime

To model the impact of staged rollout on downtime, we assume that failure in
the $Dev$ state does not incur downtime, since only internal testing is
performed at this stage. However, downtime incurred in the staged rollout and
$Ops$ states are proportional to the fraction of the user base multiplied by
the mean time to resolve defects such that the accumulated downtime increases
by $p_{i_{1}}*MTTR$ or $p_{Ops}\times MTTR$. This simplifying assumption may
be conservative, since not all users exposed to the functionality will
necessarily experience the failure. Similar to delivery time, the downtime
experienced must be modeled from staged rollout data and has important
implications for identifying an optimal deployment policy for transition
times, since conservative assumptions may unnecessarily increase delivery
time. Thus, our preliminary model expresses the total downtime as the weighted
sum $MTTR\times\sum_{i=1}^{n}p_{s(i)}$, where $s(i)$ denotes the state in
which the $i^{th}$ failure occurs.

####  2.2.3. Reward Modeling

To model staged rollout as a reinforcement learning problem, this section
defines the reward function according to techniques from multi-objective
reinforcement learning (Liu et al., 2014). Specifically, a linear combination
of the delivery and down time measures defined in Sections 2.2.1 and 2.2.2 at
each time step. To facilitate the selection of weights in an intuitive manner,
our ongoing research seeks to specify constrained optimization problems to
natural forms such as (i) minimize delivery time while limiting downtime to a
specified constraint or (ii) minimizing downtime while limiting delivery time
to a specified constraint.

##  3\. Algorithms

This section describes the Q-learning (Watkins, 1989) algorithm to estimate
the rewards of state-action pairs, the upper confidence bound exploration
strategy, and a naive policy enumeration approach to establish a baseline with
which to compare the performance of reinforcement learning.

###  3.1. Q-Learning

Reinforcement learning (Sutton and Barto, 2018) employs a Markov decision
process to estimate the reward of each state-action pair denoted $E[r(s,a)]$,
which determines a policy $\pi(s)$ that informs the best action to take in any
given state. Q-learning (Watkins, 1989) is a model free algorithm to estimate
the reward of a state-action pair according to the following update equation

(1)
$Q(s_{t},a)=Q(s_{t},a)+\alpha(r+\gamma\max_{a^{\prime}\in\mathbb{A}}Q(s_{t+1},a^{\prime})-Q(s_{t},a))$

where $s_{t}$ is the state at time step $t$, $a$ the action taken,
$\alpha\in(0,1)$ the learning rate, $r$ the reward, $\gamma\in(0,1)$ the time
discount, and $a^{\prime}$ the action in the next state ($t+1$) chosen from
all available actions ($\mathbb{A}$). A time discount close to one expresses a
willingness to defer rewards farther into the future, whereas small values of
gamma increase preference for immediate rewards. The learning rate $\alpha$ is
a multiplier which determines how much emphasis to place on the most recently
observed reward when updating the reward estimate of a state-action pair.

###  3.2. Upper Confidence Bound (UCB)

The upper confidence bound exploration strategy combines the $Q$-value with an
adaptation from statistical upper confidence limits to choose the next action
in a given state.

(2)
$a^{\prime}=\max_{a\in\mathbb{A}}Q(s,a)+c\times\sqrt{\log\left(\frac{t+1}{n(s,a)}\right)}$

where $c>0$ is a scaling constant, $t$ is the present time step, and $n$ is
the number of times that a state-action pair has been selected previously.
Thus, the UCB exploration strategy selects an action according to the
estimated reward as well as the uncertainty embodied in the finite sample size
associated with state-action pairs. The upper confidence term is larger for
smaller values of $n$, encouraging exploration of less frequently taken state-
action pairs. Small values of $c$ tend toward exploitation of the action with
the highest $Q$-value, whereas large values of $c$ encourage exploration of
the least frequently taken actions.

###  3.3. Naive Policy Enumeration

Naive policy enumeration does not use reinforcement learning, is not capable
of being used in an online fashion, and is therefore not a competitive
alternative. It is a deterministic approach based on emulation with the full
data to approximate the Pareto optimal front of tradeoffs between delivery
time and downtime from a range of stationary policies in order to objectively
compare the performance of a reinforcement learning approach in terms of its
distance from optimality. Specifically, state transition policy vector
$\bm{\lambda}=\langle\lambda_{Dev,i_{1}},\lambda_{i_{1},Ops}\rangle$ for
Figure 1 is applied with a range of deterministic values such as the cross
product of policies, where $\lambda_{Dev,i_{1}}$ and $\lambda_{i_{1},Ops}$
respectively denote the time to spend in the $Dev$ and staged rollout state
without failure before transitioning to the staged rollout or $Ops$ state.
Thus, naive enumeration applies a range of deterministic policy vectors on the
same data set and computes a tuple including the downtime and delivery time
for each policy vector. A plot consisting of only the tuples that correspond
to non-dominated policies (lowest delivery time for a fixed downtime or lowest
downtime for a fixed delivery time) constitute the Pareto optimal front for
comparison with reinforcement learning strategies.

##  4\. Metrics

To impose greater rigor and promote objective comparison of reinforcement
learning approaches, we define two quantitative metrics, including range,
which measures the flexibility of an approach, and average suboptimality,
which may be regarded as a performance-oriented metric. Both of these metrics
are computed relative to a baseline Pareto optimal curve determined from naive
policy enumeration. The discussion that follows introduces these metrics,
provides a mathematical formula, and explains what it quantifies as well as
how the metric supports objective comparison.

###  4.1. Range

The range of an approach ($a$) with respect to an objective ($o$) is

(3) $range_{a,o}=\frac{\max{a,o}-\min{a,o}}{\max{naive,o}-\min{naive,o}}.$

Thus, Equation (3) is the range of values an approach achieves with respect to
an objective divided by the range of the naive approach, which is treated as a
baseline for the sake of normalization. The range of the naive approach is
therefore $1.0$. Approaches that attain a range greater than $1.0$ for an
objective produce a wider Pareto front than the naive approach, which directly
impacts policy selection, since it will not be possible to identify a Pareto
optimal policy below the minimum or above the maximum of a specific objective
with an approach.

###  4.2. Suboptimality

The suboptimality of a specific policy obtained by an approach with respect to
an objective is

(4) $subopt_{p,a,o}=\frac{o_{p,a}}{o_{p,naive}}$

where $o_{p,a}$ is the value of the objective achieved by policy $p$ obtained
with approach $a$ and $o_{p,naive}$ is the value of the objective achieved by
a policy determined from the naive approach. Thus, Equation (4) is simply the
ratio of the value achieved by an approach for a specified objective divided
by an equivalent value for the naive approach when all other objectives are
held constant. In other words, we divide the value associated with a point on
the Pareto curve of an approach by an equivalent point on the curve for the
naive approach. Often, there is no naive policy that possesses the exact same
value as the point for the policy under consideration. We therefore
approximate suboptimality by linearly interpolating two points on the naive
Pareto curve in order to obtain a precisely matching value. The average
suboptimality ($E[subopt_{p,a,o}]$) is a summary statistic, which computes the
supoptimality for each point (individual policy) of a specified approach and
objective and then calculates the average of these values.

##  5\. Results

To demonstrate the potential to automate staged rollout with reinforcement
learning, the SYS1 data set (Lyu, 1996), which documents $n=136$ times at
which unique defects were detected during approximately 25 hours of testing
was employed. The weight on delivery time ($w_{0}$) was varied in the interval
$(0,1)$, while the remaining weight was placed on downtime such that
$w_{1}=1-w_{0}$. For each combination of weights, Q-learning with the upper
confidence bound exploration strategy given in Equation (2) was applied with
parameters $c=0.15$, learning rate $\alpha=0.15$, and time discount
$\gamma=0.999999$. It should be noted that Q-learning with the UCB strategy
learned as progress was made along the testing timeline of the SYS1 data set.
Therefore, this method did not require separate data sets for training. Naive
policy enumeration was performed with pairs of values from the cross product
of $t_{Dev,i_{1}}=\\{1,100,200,\dots,10000\\}$ and
$t_{i_{1},Ops}=\\{1,100,200,\dots,10000\\}$. For the sake of illustration,
parameters of the staged rollout model were set to $MTTR=10$, $n_{Dev}=50$,
$n_{i}=1000$, and $n_{Ops}=10000$.

Figure 2 compares the range of policies (points) identified by UCB as well as
those determined by naive policy enumeration, which can only serve as a
baseline because it is deterministic, unlike reinforcement learning, which
makes decisions as data becomes available. Figure 2 indicates that UCB
policies with low delivery times and downtime greater than $300$ are
competitive with naive enumeration. However, there is a visible gap between
the downtime achieved by UCB compared to naive enumeration for downtimes less
than $300$. One possible explanation is that UCB must perform some amount of
exploration. While parameter tuning could certainly reduce this gap, a more
important enhancement will be to incorporate software engineering specific
testing factors that reinforcement learning can use to drive staged rollout
decisions. Moreover, one may regard the results obtained by naive policy
enumeration as “lucky”, since they simply specify fixed values for times to
wait before transitioning between states when no failure is experienced. The
average suboptimality metric given in Equation (4) summarizes this gap as a
single number. For example, the average suboptimality of UCB with respect to
downtime and delivery time were $2.79$ and $2.72$, indicating the potential
margin for improvement.

Figure 2. Tradeoff between downtime and delivery time on SYS1 data

[]

Figure 2 also shows that the best policy to minimize downtime that UCB could
identify was approximately $70$ (upper left most grey dot on Pareto curve).
Thus, naive enumeration was able to achieve less downtime at the expense of
higher delivery times. The range metric given in Equation (3) summarizes the
width of the UCB policies relative to naive enumeration, which were $0.80$ and
$0.83$ for downtime and delivery time respectively, indicating that policies
identified by UCB were approximately $80\%$ as wide as the corresponding naive
policies.

##  6\. Future Plans

Our preliminary results obtained for Q-learning combined with the upper
confidence bound approach suggest that reinforcement learning is a viable
approach to automate staged rollout. Our ongoing work will consider
alternative state of the art reinforcement learning methods to match or exceed
naive policies and provide a greater range of flexibility. Techniques to
specify a desired balance of objectives as a constrained optimization problem
will provide a layer of abstraction to automatically select and apply a
suitable policy on the Pareto front. Toward this end, simple and efficient
techniques such as binary search can be used to identify a policy that best
matches stakeholder needs.

## Acknowledgment

This material is based upon work supported by the National Science Foundation
under Grant Number 1749635 and the Office of Naval Research under Award Number
N00014-22-1-2012. Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the authors and do not necessarily
reflect the views of the National Science Foundation or Office of Naval
Research.

## References

  * (1)
  * Abdallah and Kaisers (2016) Sherief Abdallah and Michael Kaisers. 2016.  Addressing environment non-stationarity by repeating Q-learning updates.  _The Journal of Machine Learning Research_ 17, 1 (2016), 1582–1612.
  * Brunnert et al. (2015) Andreas Brunnert, André van Hoorn, Felix Willnecker, Alexandru Danciu, Wilhelm Hasselbring, Christoph Heger, Nikolas Herbst, Pooyan Jamshidi, Reiner Jung, Joakim von Kistowski, et al. 2015\.  Performance-oriented DevOps: A research agenda.  _arXiv preprint arXiv:1508.04752_ (2015).
  * Farr (1996) W. Farr. 1996.  _Handbook Of Software Reliability Engineering_.  McGraw-Hill, New York, NY, chapter Software Reliability Modeling Survey, 71–117.
  * Khetarpal et al. (2020) Khimya Khetarpal, Matthew Riemer, Irina Rish, and Doina Precup. 2020\.  Towards continual reinforcement learning: A review and perspectives.  _arXiv preprint arXiv:2012.13490_ (2020).
  * Krass (1990) Dmitry Krass. 1990\.  _Contributions to the theory and applications of Markov decision processes_.  Ph.D. Dissertation. The Johns Hopkins University.
  * Liu et al. (2014) Chunming Liu, Xin Xu, and Dewen Hu. 2014.  Multiobjective reinforcement learning: A comprehensive overview.  _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ 45, 3 (2014), 385–398.
  * Lo and Huang (2006) J. Lo and C. Huang. 2006\.  An integration of fault detection and correction processes in software reliability analysis.  _Journal of Systems and Software_ 79, 9 (2006), 1312–1323.
  * Lyu (1996) Michael Lyu. 1996\.  _Handbook of software reliability engineering_ (2 ed.).  IEEE computer society press CA.
  * Nafreen et al. (2020) Maskura Nafreen, Melanie Luperon, Lance Fiondella, Vidhyashree Nagaraju, Ying Shi, and Thierry Wandji. 2020.  Connecting Software Reliability Growth Models to Software Defect Tracking. In _International Symposium on Software Reliability Engineering_. IEEE, 138–147.
  * Nelson (2009) Wayne Nelson. 2009\.  _Accelerated testing: statistical models, test plans, and data analysis_. Vol. 344.  John Wiley & Sons.
  * Sutton and Barto (2018) Richard S. Sutton and Andrew G. Barto. 2018.  _Reinforcement Learning: An Introduction_.  A Bradford Book, Cambridge, MA, USA.
  * Tarvo et al. (2015) Alexander Tarvo, Peter F Sweeney, Nick Mitchell, VT Rajan, Matthew Arnold, and Ioana Baldini. 2015\.  CanaryAdvisor: a statistical-based tool for canary testing. In _International Symposium on Software Testing and Analysis_. 418–422.
  * Vamplew et al. (2017) Peter Vamplew, Rustam Issabekov, Richard Dazeley, Cameron Foale, Adam Berry, Tim Moore, and Douglas Creighton. 2017.  Steering approaches to Pareto-optimal multiobjective reinforcement learning.  _Neurocomputing_ 263 (2017), 26–38.
  * Van Moffaert et al. (2013) Kristof Van Moffaert, Madalina M Drugan, and Ann Nowé. 2013\.  Hypervolume-based multi-objective reinforcement learning. In _International Conference on Evolutionary Multi-Criterion Optimization_. Springer, 352–366.
  * Velayutham (2021) Girish Velayutham. 2021\.  Artificial Intelligence assisted Canary Testing of Cloud Native RAN in a mobile telecom system.
  * Watkins (1989) Christopher John Cornish Hellaby Watkins. 1989\.  Learning from delayed rewards.  (1989).
  * Zhao et al. (2018) Zhenyu Zhao, Mandie Liu, and Anirban Deb. 2018.  Safely and quickly deploying new features with a staged rollout framework using sequential test and adaptive experimental design. In _International Conference on Computational Intelligence and Applications_. IEEE, 59–70.


# Grid-based state space exploration for molecular binding

Hana Zupan  Frederick Heinz  Bettina G. Keller<EMAIL_ADDRESS>Department of Biology, Chemistry, Pharmacy, Freie Universität Berlin,
Arnimallee 22, 14195 Berlin, Germany

###### Abstract

Binding processes are difficult to sample with molecular-dynamics (MD)
simulations. In particular, the state space exploration is often incomplete.
Evaluating the molecular interaction energy on a grid circumvents this problem
but is heavily limited by state space dimensionality. Here, we make the first
steps towards a low-dimensional grid-based model of molecular binding. We
discretise the state space of relative positions and orientations of the two
molecules under the rigid body assumption. The corresponding program is
published as the Python package molgri. For the rotational component of the
grids, we test algorithms based on Euler angles, polyhedra and quaternions, of
which the polyhedra-based are the most uniform. The program outputs a sequence
of molecular structures that can be easily processed by standard MD programs
to calculate grid point energies. We demonstrate the grid-based approach on
two molecular systems: a water dimer and a coiled-coil protein interacting
with a chloride anion. For the second system we relax the rigid-body
assumption and improve the accuracy of the grid point energies by an energy
minimisation. In both cases, oriented bonding patterns and energies confirm
expectations from chemical intuition and MD simulations. We also demonstrate
how analysis of energy contributions on a grid can be performed and
demonstrate that electrostatically-driven association is sufficiently resolved
by point-energy calculations. Overall, grid-based models of molecular binding
are potentially a powerful complement to molecular sampling approaches, and we
see the potential to expand the method to quantum chemistry and flexible
docking applications.

molecular association, state space discretisation, rotation grids, scientific
programming, coiled-coils

††preprint: AIP/123-QED

##  I Introduction

Molecular binding, i.e. the process in which two or more molecules form a
long-lived complex, is ubiquitous but very difficult to study, both
experimentally and computationally. Binding processes range across several
length and time scales. Examples include dimers of small solvent molecules
such as H2O dimers or HF dimers [1], solvation shells around ions in solution
[2], chelate-complexes that can mask more complex organic molecules [3],
association of molecules on two-dimensional surfaces including biological
membranes, binding of small molecular inhibitors to protein receptors [4, 5],
the formation of protein-protein complex [6, 7], and self-assembly of
molecular monomers into almost macroscopic structure such as hydrogels or
fibrils [8, 9, 10].

Using molecular-dynamics (MD) simulations combined with rare-event techniques
[11, 12, 13, 14], one can now model and analyze molecular binding at atomistic
detail [15, 16, 15]. However, MD simulations rely on sampling the
conformational space and assume that the ergodic hypothesis is fulfilled. This
is problematic for binding processes, because binding often occurs via
multiple distinct binding pathways. These pathways predominantly manifest when
the molecules are already in direct contact [4, 17, 18, 19]. But even with
large distances between associating molecules their long-range electrostatic
interactions are not isotropic but can “guide” them towards forming a complex
[20, 21, 6]. Exhaustively sampling all relative positions and relative
orientations of the binding partners is challenging, and even determining
whether or not the simulation has converged is difficult [22].

An alternative to sampling is to discretize the conformational space on a grid
and to explore its details by changing the resolution of the grid. Recently,
methods that use this approach to calculate grid-based models of
conformational dynamics have been proposed [23, 24, 25, 26, 27]. Selecting the
grid points so that interactions in all relevant degrees of freedom are
resolved uniformly is a crucial step in this process. We here make the first
steps towards grids for the binding process of two molecules. As a first
approximation, we model each of the molecules as a rigid body. Fixing the
first molecule at the origin of the coordinate system, the relative positions
and orientations of the second molecule can be described by a combination of
three subspaces: ($i$) translation in radial direction, ($ii$) rotation around
the origin, and ($iii$) rotation around its centre of mass (Fig. 1.A).
Discretising each of the three subspaces yields a regular grid which we call a
trans-rot-rot grid. Each of the grid-points can then be transformed into a set
of matrices acting on position vectors of the atoms that yield a specific
relative arrangement of the two molecules, i.e. a molecular structure.

A critical point in this process is the construction of uniform rotational
grids. Since rotations can be represented as points on a unit sphere (Fig.
1.B), this problem is mathematically closely intertwined with uniformly
distributing $N$ points on a sphere. We test rotation grids based on random
and systematic Euler angles [28], random quaternions [29], polyhedra [30, 31]
and polytopes [32].

For the construction and analysis of rotational grids, we can build on a vast
body of previous research [33] since most areas of scientific computing need
to consider grids on spheres or rotation spaces (e.g. geo[34] and atmospheric
[35] sciences, robotics [36], computer animations [29]). Rotation matrices and
the corresponding quaternion representations of rotations are also frequently
used in the context of molecular sampling, most notably for the superposition
of molecular structures [37]. Other applications include integrators of the
equations of motion for rigid bodies [38, 39, 40], reverse mapping in
multiscale simulations [41], analysis of rotational distributions [42],
rotational entropies [43] or rotational correlation functions [44].

We implemented trans-rot-rot discretisation of association space in the
program molgri and provide it as a free Python package. The program takes
structure files of two molecules as inputs and returns the relative positions
and orientations of the molecules on a trans-rot-rot grid. The output is
written in the format of a MD-trajectory, so that the grid point energies can
be calculated using a standard MD program.

In this contribution, we briefly discuss the theory of rotations and explain
the implementation and the usage of the program molgri. Next, we test
rotational grids generated by six different algorithms regarding their
uniformity and suitability for molecular studies. We then demonstrate how
trans-rot-rot grids can be used to analyze a molecular binding process. We
discuss the position- and orientation-dependent contributions to the
interaction energies of two water molecules, i.e. molecules which are modelled
as rigid bodies in classical force fields. As a second example, we consider
the interaction between a Cl--anion and a coiled-coil peptide [10].

##  II Theory

We introduce three common descriptions of rotations and briefly review how
they relate to each other. We also discuss measures of uniformity applicable
to rotation grids. For a systematic overview of rotation representations and
conversions between different formats we recommend a comprehensive overview
article by Diebel [33].

###  II.1 Rotation matrices

A 3D rotation matrix $R$ is an orthogonal $3\times 3$ matrix, i.e. a matrix
which fulfills $RR^{\mathrm{T}}=R^{\mathrm{T}}R=I$, where $R^{\mathrm{T}}$ is
the transpose of $R$ and $I$ is the $3\times 3$ identity matrix. When applied
to two vectors $\mathbf{v}^{\prime}=R\mathbf{v}$ and
$\mathbf{u}^{\prime}=R\mathbf{u}$, an orthogonal matrix preserves the length
of these vectors and the angle between them, which are two key properties of
rotations. The determinant of an orthogonal matrix always equals $\pm 1$. We
define orthogonal matrices with $\det R=1$ as those describing proper
rotations while those with $\det R=-1$ describe improper rotations or
rotation-inversions.

###  II.2 Euler angles

Euler angles, first introduced by Leonhard Euler [28], are based on the idea
that ($i$) any two Cartesian axes in 3D space span a plane in which elementary
rotations can be parametrized with a single angle, and ($ii$) any rotation in
3D space can be written as a sequence of three such elementary rotations
around non-repeating axes.

The triple of angles $\mathbf{u}_{xyz}=[\phi,\theta,\psi]^{T}$ is called an
Euler angle set and is a description of a 3D rotation in which a vector
undergoes a sequence of three elementary (planar) rotations: first around axis
$z$ for an angle $\psi$, then around $y$ for $\theta$ and finally around $x$
for $\phi$. Other sequences of axes are also possible, another common choice
is the $zxz$ set.

Since the rotation axes for the elementary rotations are the axes of the
Cartesian coordinate system, their corresponding elementary rotation matrices
are simple trigonometric functions of the rotation angle. For rotation of
$\phi$ around an x-axis we have:

(1) $R_{x}(\phi)={\left[\begin{array}[]{ccc}1&0&0\\\
0&\cos(\phi)&\sin(\phi)\\\ 0&-\sin(\phi)&\cos(\phi)\end{array}\right]},$

and the $R_{y}$ and $R_{z}$ elementary matrices are simple permutations of the
matrix above [45].To transform an Euler angle set into a single rotation
matrix, we multiply the corresponding elementary rotation matrices in the
correct order:

(2) $R_{xyz}(\phi,\theta,\psi):=R_{x}(\phi)R_{y}(\theta)R_{z}(\psi).$

###  II.3 Quaternions

Another common representation of rotations is that of unit quaternions.
Quaternions have been invented by William Rowan Hamilton [46, 47] and are
mostly known for their application in space rotations. A quaternion
$\mathbf{q}$ is a 4-dimensional complex number

(3)
$\mathbf{q}=\left[q_{0},q_{1},q_{2},q_{3}\right]^{T}=a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k},$

where a, b, c, and d are real numbers and $\mathbf{i}$, $\mathbf{j}$, and
$\mathbf{k}$ are the basic quaternions that can be interpreted as unit vectors
along three perpendicular coordinate axes with the property
$\mathbf{i}^{2}=\mathbf{j}^{2}=\mathbf{k}^{2}=-1$. The length of a quaternion
is defined as

(4) $\|\mathbf{q}\|=\sqrt{q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2}}.$

A unit quaternion is a quaternion with $\|\mathbf{q}\|=1$. We state, but omit
the proof, that there always exist exactly two unit quaternions that map to
one rotation matrix. For the derivation see Refs. 46, 47, 33. A rotation of
any vector $\mathbf{v}$ to $\mathbf{v}^{\prime}$ can be described in terms of
quaternion multiplication. For this, the three-dimensional vector
$\mathbf{v}=[v_{1},v_{2},v_{3}]^{T}$ is interpreted as a quaternion in which
the first element (the scalar element) is zero:
$\mathbf{v}=0+v_{1}\mathbf{i}+v_{2}\mathbf{j}+v_{3}\mathbf{k}$. The rotation
then is

(5) $\mathbf{v}^{\prime}=\mathbf{q}\mathbf{v}\mathbf{q}^{-1},$

where $\mathbf{q}$ is a unit quaternion and

(6)
$\mathbf{q}^{-1}=\frac{\overline{\mathbf{q}}}{\|\mathbf{q}\|},\qquad\qquad\overline{\mathbf{q}}=\left[q_{0},-q_{1},-q_{2},-q_{3}\right]^{T}.$

is its reciprocal. The multiplication of quaternions in eq. 5 is defined by
the Hamilton product [46, 47, 33].

Quaternions have proven to be the most elegant solution in computational
implementation of rotations. Since there is a natural and smooth mapping from
the space of quaternions to the space of rotations (the $SO(3)$ group),
interpolation of rotations can be easily performed in quaternion space and no
singularities occur within this description, as opposed to the Euler angles
parametrisation.

###  II.4 Rodrigues’ formula

Rotational matrices can be inferred from an initial vector $\mathbf{v}$ and
rotated vector $\mathbf{v}^{\prime}$ using Rodrigues’ formula that has been
independently discovered by Euler [48] as well as Rodrigues [49]. For the
derivation, see chapter 3 of Ref. 50. In essence, we can calculate the unit
vector $\bm{\hat{\omega}}$ along the axis of a rotation that transforms
$\mathbf{v}$ into $\mathbf{v}^{\prime}$ and the angle of rotation $\theta$ as

(7)
$\bm{\hat{\omega}}=\frac{\mathbf{v}\times\mathbf{v}^{\prime}}{\|\mathbf{v}\|\|\mathbf{v}^{\prime}\|},\quad\sin(\theta)=\|\bm{\hat{\omega}}\|,\quad\cos(\theta)=\frac{\mathbf{v}\cdot\mathbf{v}^{\prime}}{\|\mathbf{v}\|\|\mathbf{v}^{\prime}\|}.$

We additionally define the skew-symmetric matrix of a vector
$\bm{\omega}=\left[\begin{array}[]{lll}\omega_{1}&\omega_{2}&\omega_{3}\end{array}\right]^{\mathrm{T}}\in\mathbb{R}^{3}$
as

(8) $[\omega]_{\times}=\left[\begin{array}[]{ccc}0&-\omega_{3}&\omega_{2}\\\
\omega_{3}&0&-\omega_{1}\\\ -\omega_{2}&\omega_{1}&0\end{array}\right].$

The transformation $\mathbf{v}\rightarrow\mathbf{v}^{\prime}$ is then
performed with

(9)
$\mathbf{v}^{\prime}=e^{[\hat{\omega}]_{\times}\theta}\mathbf{v}=(I+\sin\theta[\hat{\omega}]_{\times}+(1-\cos\theta)[\hat{\omega}]_{\times}^{2})\mathbf{v},$

which is known as Rodrigues’ rotation formula. Writing the skew-symmetric
matrices in eq. 9 in matrix form, we can directly convert this representation
into a rotation matrix multiplication.

###  II.5 Grid coverage and uniformity

We use different parameterisations of rotations to create rotation grids. In
order to compare how well the generated set of grid points covers the entire
space of rotations, we introduce a measure of uniformity. Consider a rotation
grid with $N$ grid points, each represented by a point on a unit sphere
(radius $r=1$, total surface area $A=4\pi$). We calculate the local grid
density in various sectors of the unit sphere by counting the number of grid
points $N(S(\mathbf{v},\alpha))$ within a well-defined sub-area
$S(\mathbf{v},\alpha)$ of the sphere surface and calculating the coverage
ratio

(10) $\displaystyle c(\mathbf{v},\alpha)$ $\displaystyle=$
$\displaystyle\frac{N(S(\mathbf{v},\alpha))}{N}.$

The vector $\mathbf{v}$ denotes the center of $S$, and the parameter $\alpha$
determines its surface area $A_{S}$. In a fully uniform grid, the coverage
would be equal for all sub-areas of equal size $A_{S}$ and further equal to
the ratio $A_{S}/A$.

Figure 1: A. Three types of motion in molecular binding: radial translation,
relative rotation and internal rotation of one molecule around its own centre
of mass. B. Example of a rotation grid with $N_{rot}$ grid points. C.
Measuring local grid density. Spherical cap with area $S(\mathbf{v},\alpha)$
(shaded) is defined by an axis vector v and an angle $\alpha$. The ratio of
grid points within spherical cap $N(S(\mathbf{v},\alpha))$ with $N_{rot}$ is
defined to be local grid density.

A suitable choice of sub-areas is a class of curved surfaces called spherical
caps[51], defined by the intersection of the unit sphere with a plane (Figure
1C). The intersection itself is a circle, called the base of the spherical
cap. Here, $\mathbf{v}$ is the unit vector that connects the center of the
unit sphere with the center of the cap and defines the center of $S$. The
parameter $\alpha$ is determined by the opening angle of the spherical cap,
i.e. the angle $\alpha$ between $\mathbf{v}$ and a unit vector that points to
rim of the base, and is the quantity defining the size of $S$. The area of a
spherical cap is [51]

(11) $A_{S}=2\pi(1-\cos\alpha).$

And so the expected coverage for a uniform grid is

(12) $\displaystyle c_{\mathrm{uniform}}(\alpha)$ $\displaystyle=$
$\displaystyle\frac{A_{S}}{A}=\frac{1}{2}(1-\cos\alpha)$

See Refs. 31, 30 for a more detailed discussion on measures of uniformity.

##  III Software

The Python package molgri is freely available within the framework of Python
Package Index [52]. Documentation and source code are available via GitHub
[53]. Dependencies are listed in section I of the supplementary information.

The package molgri implements a range of generation, analysis and plotting
tools. A short user guide is provided in Section III.3. As the workflow in
Figure 2 indicates, the main output of molgri are the of relative positions
and orientations of two molecules on a trans-rot-rot grid. The output is
written in the file-format of a GROMACS molecular-dynamics trajectory (.gro-
file), such that is can be further processed by standard MD programs to obtain
and analyse point energies. In this contribution, we use GROMACS for the
energy calculation, but since the .gro-file format is a standard file-format
for MD trajectories and can readily be converted in a other formats, the user
is free to use the MD program of their choice for this step. We use the term
pseudo-trajectory for the output, because it has the file format of a MD
trajectory (.gro file), but it represents a trans-rot-rot grid and not the
result of MD sampling.

molgri generates pseudo-trajectories in two steps. In the first step, the
program generates a rotation grid, that is, a collection of $N_{\mathrm{rot}}$
points distributed on a unit sphere. We implemented and tested six algorithms
for this step (see section III.1). In the second step, the rotation grid is
combined with a translation grid with $N_{\mathrm{trans}}$ elements. Overall
we obtain $N_{\mathrm{trans}}\cdot N_{\mathrm{rot}}\cdot N_{\mathrm{rot}}$
relative positions and orientations for the two-molecule system (see section
III.2).

We took particular care to build a modular and flexible program structure.
Therefore, each step outputs a standard file format that can be reused by
other applications. This also means that a rotation grid of particular size
can be reused for any pair of molecules.

Rotation gridGrid algorithmNPseudotrajectoryTranslation gridMolecular
structuresMD softwareForce fieldTopology.npy file.gro file

Figure 2: Workflow of grid-based state space exploration. The first two stages
are performed with the molgri package, the resulting pseudo-trajectories can
be further evaluated by most MD programs.

###  III.1 Step 1: rotational grid algorithms

The six algorithms we implemented to generate rotational grids are: Systematic
Euler grid, Random Euler angles, Random unit quaternions [29], 3D cube grid
[30], Icosahedron grid [31], and 4D cube grid [32]. The algorithms are
detailed in section II of the supplementary information, but we sketch the
ideas behind them here.

Systematic Euler grid is based on a uniformly spaced grid between 0 and $2\pi$
for each of the three Euler angles using the $xyz$ sequence of axes. Random
Euler angles algorithm generates a random rotation by selecting a uniform
random number in the $[0,2\pi]$ interval for each of the three Euler angles.
Random unit quaternions algorithm generates random rotations by selecting
random quaternions using the subgroup algorithm [29]. Icosahedron grid and 3D
cube grid use as a regular grid on the surface of a polyhedron to approximate
a uniform grid on a unit sphere. This approach starts with inscribing a 3D
polyhedron (a cube[30] or an icosahedron[31]) into a unit sphere, then
constructing a grid on the faces of the polyhedron and projecting the grid
points onto the unit sphere. The number of grid points $N_{\mathrm{rot}}$ can
be increased by subdividing the faces of the polyhedron. We subdivide square
faces (cube) on a layered diagonal square lattice and triangle faces
(icosahedron) on a layered triangle lattice. 4D cube grid[32] extends this
concept to polytopes (generalization of polyhedra to higher dimensional
spaces) and quaternions. In this case, a four-dimensional hyper-cube is
inscribed into a four-dimensional unit hyper-sphere, a grid is constructed on
its faces, and the grid points are projected onto the hyper-sphere. The
resulting four-dimensional vectors are interpreted as quaternions.

Some of the six algorithms generate a set of Euler angles, others a set of
quaternions and still others a set of points on a sphere. We unify all outputs
in the latter format by letting the quaternions or Euler angles act on a unit
vector in the $z$-direction and saving the rotated vector as a three-
dimensional grid point (algorithm I in section II of the supplementary
information).

Systematic Euler grid and polyhedron/polytope-based grid algorithms use
regular grids to generate the rotations, and can therefore natively only
generate rotational grids with specific numbers of grid points
$N_{\mathrm{rot}}$. To provide complete flexibility in the number of grid
points $N_{\mathrm{rot}}$, we first generate a grid of size $M\geq
N_{\mathrm{rot}}$, then order all points in a way that maximizes coverage and
truncate at $N_{\mathrm{rot}}$ (algorithm II in section II of the
supplementary information).

###  III.2 Step 2: from rotation grids to pseudo-trajectories

In the second step, the translational grid, the rotational grid of relative
positions and the rotational grid of relative orientations are combined. We
use a linear equidistant grid for the translations, and currently the same
rotational grid for selection of relative positions and relative orientations.
Simplifications due to symmetry are not taken into consideration. A full grid
in state space is a product (meshgrid) of all three component grids.

To generate a pseudo-trajectory, the rotational grid is first transformed into
a series of rotation matrices using Rodrigues formula (Equation 9).
Specifically, the unit-vector in the $z$-direction is the initial vector
$\mathbf{v}$ and the grid point is the rotated vector $\mathbf{v}^{\prime}$.
Rotation axis and angle can then be calculated using Equation 7.

The pseudo-trajectory is generated by keeping molecule 1 fixed at the origin
of the three-dimensional coordinate system while molecule 2 is translated and
rotated. Initially, the two molecules are positioned in their starting
orientations at distance $r_{0}$ (first value of the translational grid). The
$N_{\mathrm{rot}}$ rotation matrices are applied to the second molecule as
rotations around the origin and the resulting configurations are recorded. In
each of the configurations, the $N_{\mathrm{rot}}$ rotation matrices are then
applied to the second molecule as rotations around the center of mass,
generating $N_{\mathrm{rot}}$ relative orientations. The process is repeated
for each radius provided by the translational grid. All resulting
configurations are recorded in a single .gro file.

###  III.3 Using the molgri package

In this section, we provide short instructions for the installation of molgri
package and its use from command line. Assuming a compatible Python
installation is already present, the package is installed with the command

⬇

pip install molgri

After installation, molgri can be imported as a package into a Python program
as import molgri. Alternatively, we also provide three scripts that can be run
directly from the command line: 1) molgri-io that creates a standard tree of
input/output directories, 2) molgri-grid that generates rotation grids
including their plots, animations and statistical analyses and 3) molgri-pt
that generates pseudotrajectory files. All three scripts can be run with an
optional flag \--help that returns a short user guide to the script.

Running the script molgri-io with the optional flag \--examples

⬇

molgri-io –examples

provides example input .gro files including all molecules and ions used in
Results section of this contribution and generates the following directory
structure

⬇

.

|-output

| |-grid_files

| |-pt_files

| |-figures

| |-animations

| |-statistics_files

|-input

| |-CL.gro

| |-H2O.gro

| |-NA.gro

| |-NH3.gro

| |-example_protein.gro

Instead of running this script, the user can also manually create an input/
directory and copy any .gro files for which pseudo-trajectories should be
generated there.

The second script molgri-grid implements step 1 from Figure 2. It is necessary
to specify the number of grid points -N and the algorithm -algo (options:
systemE, randomE, randomQ, cube4D, cube3D, ico) to generate a rotation grid.
In addition, optional flags to this script are: \--draw that saves a plot of
grid points, \--animation for the corresponding 3D animation, \--statistics
that generates a number of files and figures analysing the uniformity and
convergence of this grid, and \--readable that saves the grid points in a
human-readable .txt format in addition to the standard .npy format. The
generated files are saved to the correspondingly named subfolder within the
output folder. For example,

⬇

molgri-grid -N 250 -algo ico –draw

–animate –statistics –readable

generates a rotational grid with 250 grid points using the algorithm
Icosahedron grid, draws and animates the grid, calculates the grid statistics,
and additionally stores the grid in human-readable format.

The last command line script molgri-pt implements step 2 and creates a
pseudotrajectory. It expects two molecular structure files in the .gro-format
which should be stored in the input/ folder. Usually it is preferential to
select the smaller molecule as the -m2 option since this is the molecule
undergoing translations and rotations. In addition, the user needs to specify
the rotational grid in the form -rot algorithm_N (see algorithm names above)
and the translational grid (flag -trans). The translational grid can be
specified in one of the following formats: a list of distances,
linspace(start, stop, num) or range(start, stop, step), where the units are in
nanometers. For example, the translational grid with grid points 1.0, 1.5,
2.0, 2.5, and 3.0 nm can be specified as

  * •

-trans "(1, 1.5, 2, 2.5, 3)"

  * •

-trans "linspace(1, 3, 5)"

  * •

-trans "range(1, 3.1, 0.5)"

In summary, the command

⬇

molgri-pt -m1 NH3 -m2 H2O -rot cube3D_15

-trans ”linspace(1, 3, 5)”

generates a pseudo-trajectory of a water molecule rotating around an ammonium
molecule, where the rotational grid is generated with a cube 3D algorithm and
contains 15 points (or is read from /output/grid_files/cube3D_15.npy if
already generated). The two structure files are expected as ./input/H2O.gro
and ./input/NH3.gro, and the translational grid is the one discussed above.
Finally, the user may select an optional flag \--only_origin which suppresses
the generation of different orientations of the second molecule and only
returns $N_{\mathrm{rot}}\cdot N_{\mathrm{trans}}$ structures. This is useful
if the second molecule is spherically symmetric (e.g. a single ion or an
atom).

##  IV Results

###  IV.1 Comparison of rotational grid algorithms

Figure 3: Comparison of rotation grid algorithms. A. Rotation grids with N=600
points are shown in upper left corners. The variance in grid densities across
sphere surface is displayed as a set of violin plots. As a measure of
uniformity we select 100 random axes and measure the density of the grid in
this region by counting the number of points $N_{sc}$ that fall within an
angle $\alpha$ of the axis and calculating the coverage $N_{sc}/N$. Broad
distributions suggest that coverage is not uniform. B. The same rotation grids
are used to create pseudo-trajectories for the water-water system. Potential
energies [kJ/mol] of structures with various orientations and positions of the
second water molecule are plotted, energy minima shown in dark blue. The
Systematic Euler algorithm misses one of the three minima (below the hydrogen)
and the cube 4D algorithm covers it with only one point, demonstrating the
errors that may arise when using grids with poor uniformity.

Algorithms presented in Section III.1 provide several ways to discretise
rotation space, each generating a particular set of $N_{\mathrm{rot}}$ points
on a sphere. The insets in Figure 3 A show the grid points generated by each
of the algorithms for $N_{\mathrm{rot}}=600$. One can easily distinguish grid-
based approaches from stochastic ones by the regularity and the symmetry of
the grids. Note that in the Systematic Euler grid, the grid points align in
circles along the “latitudes” around the sphere and are denser at the
“equator” than at the “poles”. Since the random Euler angle algorithm is a
direct stochastic analogue of the Systematic Euler grid, it also over-samples
along the equator and under-samples the pole areas. Random unit quaternions
algorithm samples from uniform distribution in the unit sphere (analytical
proof in [29]), but due to fluctuations in a random sample with finite number
of points, this does not necessarily yield the most uniform grid. Polyhedra-
based algorithms (3D cube grid and Icosahedron grid) show obvious regularities
in their structure. To quantitatively compare the grids, we recall the measure
of uniformity introduced in section Section II.5 where the local density of
grid points was measured by counting the number of grid points within a
spherical cap area (Equation 11) and evaluating variations from ideal grid
density given in Equation 12. Local density is measured around 100 randomly
chosen axes and the variation of measurements shown in form of a violin plot
in Figure 3 A. The parameter $\alpha$ indicates the size of the spherical cap
area in which counting is performed. The narrower the distribution, the more
uniform the grid. The two algorithms that are based on inscribing polyhedra
into three-dimensional sphere, 3D cube grid and Icosahedron grid, produce by
far the most uniform grids across all values of values of $\alpha$.

It is important to consider several different angles to test uniformity at all
scales, from narrow areas surrounding the axis to ones covering almost entire
sphere (we test angles ranging from $\pi/6$ to $5\pi/6$). How variance changes
with $\alpha$ provides additional information. We observe, for example, that
variance in random Euler angle grid changes with $\alpha$, suggesting that the
algorithm fails to achieve true randomness (random quaternion algorithm
performs better in this regard). Using a single alpha value can also be
deceiving as it may reflect inherent symmetries in grid construction. Notable
examples are cubic grids with essentially nonexistent variance in grid
densities for $\alpha=\pi/2$. For the 3D and 4D cubic grids, the subdivision
is performed in the same way for each of the eight faces, meaning that each
spherical cap covering a quarter of the sphere ($\alpha=\pi/2$) always
contains the same collection of points ($\pm 1$ per face for $N$ not divisible
by 8).

Grid uniformity is not only of theoretical interest but has real consequences
when using rotation grids to discretise the space of molecular rotations.
Figure 3 B displays a simple molecular example where the grids discussed above
were employed to study a system of two water molecules. In this figure, areas
where structures with low potential energy were found are shown in dark blue
color. In most cases, we observe three minima: a broad one near the free
electron pairs of the central oxygen atom and two more pointed ones near each
of the central hydrogens, indicating three possible hydrogen bonding patterns.
The difference between grid algorithms becomes apparent when we note that one
of the three minima (near bottom hydrogen) is completely missing in the
Systematic Euler grid and is only represented by a single point in the 4D cube
grid. Although a very dense (600 points) grid was used, grids with poor
uniformity failed to identify all regions of interest even in this simple
example. This underlines the necessity of examining topology when selecting
parametrisation of non-Euclidean spaces.

Figure 4: Time needed to construct rotation grids using six different
algorithms. Polyhedra grids (icosahedron, cube 3D) are orders of magnitude
more time intensive than random grids (random Euler angles, quaternions), but
total time needed to create grids with up to 1000 points is still under 100
seconds.

Figure 4 compares how time-intensive the six rotational grid algorithms are.
The timings have been measured on a Debian GNU/Linux 11 operating system,
equipped with an Intel(R) Core(TM) i5-8500 CPU @ 3.00GHz and 30 GB RAM. Both
random algorithms, Random Euler angles and Random unit quaternions, are very
efficient, while the 3D cube grid and Icosahedron grid show much poorer
scaling with the number of grid points. The almost constant timing for long
intervals of $N_{\mathrm{rot}}$ originate from the fact that polyhedra-based
algorithms do not natively yield grids with arbitrary number of grid points -
larger grids are created and truncated according to an unification algorithm
found in Section II of the Supplementary information. However, even for the
more time-consuming algorithms the computational costs are in the range of
seconds for grids with up to 1000 points, which is very small compared to the
cost of the entire analysis, which is dominated by the energy calculation
along the pseudo-trajectory.

In summary, polyhedra-based grids are remarkably uniform. The somewhat larger
computational cost for generating these grids are well worth the improved
quality.

###  IV.2 Water-water system

In the second step of the molgri procedure, rotation and translation grids are
combined to systematically parameterise all three movements with which we
parameterise association space - translation, relative rotation and internal
rotation, see Figure 1 A. After selecting two molecules, a pseudo-trajectory
based on the meshgrid of the three grids is generated.

We first demonstrate the use of pseudo-trajectories on a system of two water
molecules, since this is a small, well-studied example that nonetheless
displays several oriented hydrogen-bonding patterns, enabling us to study how
radial and rotational degrees of freedom contribute to potential energy.
Moreover, rigid models of water are common in classical MD, meaning that a
direct comparison with established methods is possible.

|  Potential [$\mathrm{kJ}/\mathrm{mol}$] | H-bond length [nm]
---|---|---
systematic Euler a. | –27.67 | 0.28
icosahedron grid | –27.39 | 0.28
random Euler a. | –27.35 | 0.28
random quaternions | –27.24 | 0.27
3D cube grid | –27.23 | 0.27
MD run | –27.09 | 0.28
4D cube grid | –27.03 | 0.27
Table 1: Comparing minima of (pseudo)trajectories, each of the methods
evaluating 60.000 points in state space.

In Table 1 we perform such a comparison of minimal potential energies found in
a relatively sparse ($N_{rot}=100$, $N_{trans}=5$, 0.3-0.32 nm) pseudo-
trajectory of two rigid (TIP3P) water molecules with a standard GROMACS
trajectory with the same number of frames. Since all found minima featured a
hydrogen bond, we also compare hydrogen bond lengths. While there is some
difference between different rotation algorithms, all minima fall within 2.5%
of the comparison calculation and all but one molgri algorithms even find a
deeper minimum than a simulation with the same number of points. The plots in
previous section (Figure 3) also suggest that the minima found by molgri occur
at locations expected by partial charge distribution in water molecules.

Figure 5: A. Coulomb and LJ energies [kJ/mol] for the full pseudotrajectory
(icosahedron grid, $N_{\mathrm{rot}}=600$, $N_{\mathrm{trans}}=5$) separated
into three plots based on the the radial COM-COM distance of both water
molecules. B. Pseudo-trajectory structures with lowest (blue box) and highest
(red box) Coulomb/LJ contributions.

In the the next step, we examine how different components of potential energy,
specifically the Coulomb and Lennard-Jones contributions, vary with distance
and orientation. The ability to systematically examine spatial dependence of
energy contributions is a unique feature of our approach since classical
trajectories rarely enable us to examine the full range of states, strongly
biasing the sample towards a few easily accessible minima. For this
examination, we construct longer pseudo-trajectories of the same system using
Icosahedron grid with $N_{rot}=600$ and $N_{trans}=5$ between 0.3 and 0.32 nm,
generating $600\cdot 600\cdot 5=1.8$ million time steps in total.

In order to separate translational variation from rotational, we plot three
shells with different radial distances in separate subplots of Figure 5 A. We
plot a point at the center of mass (COM) of the second water molecule and
color-code it according to its Coulomb or Lennard-Jones energy. Since only a
subset of the 1.8 million structures can be reasonably visualised, we first
perform two levels of selection: 1) among all structures with the same COM but
different internal orientations, we select the one with lowest energy and 2)
of the points selected in step 1), 100 lowest and 100 highest energies are
plotted. On the bottom of the same Figure, 5 B, we plot the structures with
the lowest and highest energy contribution so that internal orientation of the
second water molecule can also be seen.

Completely different patterns arise for the two energy components. In the case
of Lennard-Jones potential, radial distance plays a major role (as can be
expected from the functional form of this potential) while the hydrogen-
bonding-specific patterns are not present. Hydrogen-hydrogen interaction is
the most favourable LJ interaction in this system while oxygen-oxygen one,
especially at smaller radii, is most repulsive.

On the other hand, Coulomb energies barely vary with radial distance but a
strong and persisting orientational pattern is observed instead. This could
again be predicted by the form of Coulomb potential which scales with $1/r$,
$r$ being the distance between charges. The lowest- and highest energy
structures confirm that the cause is electrostatic repel between positively
charged hydrogens and attractive interaction of a hydrogen with free electron
pairs on an oxygen.

While the contribution analysis presented here is anything but novel for this
simple system, we expect this kind of systematic study to be valuable for
complex systems featuring a mix of interaction potentials, something that is
very typical for protein systems.

###  IV.3 Polypeptide-ion system

Figure 6: A. An example of a pseudo-trajectory with a sparse grid
($N_{\mathrm{rot}}=30$, $N_{\mathrm{trans}}=8$, distances 0.5 - 2.25 nm).
Molecular structure of the coiled-coil peptide is shown in grey and the probe
Cl- ion in orange. All frames of the pseudo-trajectory are shown overlapped,
each contains one Cl- ion. In B and C parts of the figure,
$N_{\mathrm{rot}}=1000$ is used, otherwise, the set-up is equivalent. B.
Change in the potential energy of the coiled-coil-Cl- system through the usage
of energy minimisation of each frame. Potential energy before (black) and
after (red) energy minimization is shown for frames ordered from highest to
lowest energy. Additionally, the relative change in energy by minimization is
shown in blue. Some structures (zone I) exhibit immense relaxation during the
minimisation; in these cases the initial structure featured overlapping
Lenard-Jones potentials. Other structures (zone III) remain essentially
unchanged. C. Example pseudo-trajectory frames from zones I, II and III. The
initial pseudo-trajectory structure is shown in grey and the structure after
minimisation in red. Only the movement of atoms within 0.6 nm of the Cl- ion
are visualized in red; changes outside this sphere were minimal fluctuations.

The water-water dimer is an often modeled system without any internal degrees
of freedom and can be treated accurately as a system of two rigid bodies. To
test the implications of the rigid-body assumption on a system with many
internal degrees of freedom, we use molecularRotationalGrids to evaluate the
system of a chlorine anion interacting with a coiled-coil peptide (Figure 6).

The coiled-coil peptide hFF03 [10] is a parallel homo-dimer of two peptide
strands, each of which forms an $\alpha$-helix. The resulting coiled-coil
dimer is stabilized via a leucine zipper [54].The ends of the peptide strand
are not capped and the peptide is represented at pH= 7. Consequently, the amid
group of the N-terminal end of the coiled-coil dimer and amine group of the
lysine side chains are positively charged. Whereas the carboxyl group of the
C-terminal end and the glutamatic acid side chains are negatively charged.

Each peptide strand consists of 26 amino acids, and the $\alpha$-helix
measures approximately 3.8 nm in length. In the coiled-coil structures, the
helices are approximately 1.2 nm apart. In total, the system has just under
1000 atoms.

To calculate the rotational grid for coiled-coil peptide and the chlorine
anion, we placed the center of mass of the coiled-coil peptide at the origin
of the coordinate system and aligned the $\alpha$-helices with the $z$-axis.
The relative positions of the chlorine anion Cl- were generated using the
package molgri with a translational grid of 0.50, 0.75, 1.00, 1.25, 1.50,
1.75, 2.00, and 2.25 nm, i.e. an equidistant grid with $N_{\mathrm{trans}}=8$.
For the rotational grid of the relative positions, we used an Icosahedron grid
with $N_{\mathrm{rot}}=1000$ grid points. Since Cl- is a single ion, no
rotational grid for the relative orientations was needed. The combined grid
had 8000 grid points and the resulting pseudo-trajectory as many frames. To
give an impression of such a grid, Figure 6 A shows the relative positions of
a Cl- anion on a similar but more sparse grid with $N_{\mathrm{trans}}=8$ and
$N_{\mathrm{rot}}=30$.

We next calculated the interaction energy between Cl- and the coiled-coil
peptide using GROMACS rerun on the pseudo-trajectory. Since no dynamics are
present, we cannot estimate kinetic energy; the interaction energy that we
obtain consists of the Coulomb and Lennard-Jones interactions. Since the
coiled-coil peptide is rod-like and the set of Cl- configurations spherical,
some of the structures at smallest radii inevitably feature the chlorine anion
overlapping with the peptide, leading to enormous Lennard-Jones interaction
energies. We assume that many of these frames can be discarded as unphysical
cases.

To obtain more realistic grid point energies, the internal degrees of freedom
of the coiled-coil peptide can be relaxed to the nearest local minimum. We
thus ran a steepest descent energy minimization using GROMACS on every frame
of the pseudo-trajectory while constraining the position of the Cl- anion and
the peptide backbone atoms to make sure that the system does not relax into
the neighboring grid point. In three of the 8000 grid points, the energy
minimization did not lower the energy, despite very high initial energies. We
excluded these grid points from the analysis. The energy decrease for the
other grid points is shown in Figure 6 B as absolute (black and red line) and
relative (blue line) energy difference. On average, the energy after
minimisation is (2200$\pm$ 480) kJ/mol per grid point.

The changes in energy after minimisation can be roughly divided into three
categories, labeled I, II and III in Figure 6 B and C. The largest energy
differences (I) occur for grid points at which the Cl- ion overlaps with
several atoms of the coiled-coil peptide, usually because it is placed inside
the coiled-coil. Moderate energy changes (II) occur when the Cl- overlaps with
a single side-chain. In the example in the Figure 6 C-II Cl- is initially
placed on top of a lysine side chain which is located at the exterior of the
coiled-coil and can easily rotate away during energy minimization. In roughly
half of the grid points, Cl- is placed outside of the coiled-coil peptide and
only minimal adjustments are made during the energy minimization (III).

We draw the conclusion that minimization is crucial in the areas where
Lennard-Jones interaction, especially its repulsive part, is the dominant
contribution to energy. Structures from region I and II benefit from
subsequent optimization while its effect in zone III is minimal. We can
imagine that in docking studies, for example, molgri approach can be used on
its own for electrostatically-lead process of the ligand approaching the
protein while the subsequent process of the ligand docking within the protein
structure can be first studied on a grid and then followed by short
optimisation runs.

Since we are using a single anion as our association particle, we can regard
it as an electric field probe and expect it to identify positively charged
regions of state space. Figure 7 shows 100 frames of coiled-coil-Cl- pseudo-
trajectory that had lowest potential energies (before and after minimisation).
We indeed observe that electrostatic interaction is the driving force of
association since all low energy structures feature anions congregating around
the positively charged N-terminus or the amino groups of the lysine side
chains. In addition, changes to the lowest-energy structures with minimization
are small. This supports our hypothesis that molgri on its own is a suitable
tool to identify regions of interest (those featuring an appropriate charge
and a free approach vector) in association studies while energy minimization
can be used to fine tune the results.

Figure 7: 100 frames of coiled-coil-Cl- pseudo-trajectory
($N_{\mathrm{rot}}=1000$, $N_{\mathrm{trans}}=8$, distances 0.5 - 2.25 nm)
with lowest potential energies before (left) and after energy minimization
(right). The protein structure (grey) is shown in the initial (non-minimized)
form in both figures, Cl- anions are shown in orange. The red dotted spheres
represent negatively charged carboxyl groups and the blue ones represent
positively charged amino groups.

##  V Computational Methods

To test the pseudo-trajectories we generated with molgri we used GROMACS 2022
[55, 56, 57, 58] rerun function to perform point energy calculations of
configurations proposed by our algorithm. We also used GROMACS to perform
comparison MD simulations. Plots of molecular structures were obtained with
VMD 1.9 [59]. Some figures were combined and labelled in Inkscape 1.0 [60].

For comparison calculations, we used GROMACS simulations with geometry
optimisation (without the rerun option). For the water-water example system,
we used the rigid TIP3P [61] model of water as a force field and performed
simulations with 10 fs time steps. In the polypeptide example systems
(hFF03-Cl-) we used the Amber ff99SB-ILDN protein force field [62], dielectric
constant of water 78.4 (at 300K). All simulations were performed at 300 K with
particle mesh Ewald (PME) electrostatics and various trajectory lengths
(usually up to 105 time steps). All parameters were kept constant between the
pseudo-trajectory calculation and comparison MD run, although some parameters
were irrelevant for a rerun (e.g. time step).

GROMACS was also used for the frame-wise minimisation through the steepest
descent algorithm. During this process, the backbone and ion where restrained
with GROMACS restraint function, which applies a symmetric force on the chosen
atoms. This approach allows changes to the local peptide structure without
deforming the whole structure or landing in a neighbouring grid point. To
encourage the adaptation of the backbone over ion movement we applied position
restraints of 1.000.000 kJ/mol nm² to Cl- and 1000 kJ/mol nm² to protein
backbone atoms. Due to program limitations, the pseudo-trajectory needs to be
cut apart, each frame minimised individually and assembled together again.

##  VI Conclusion

We offer a fresh perspective on state space exploration for molecular binding
processes by replacing force-field-based sampling of the two-body interaction
with a regular grid. We developed a python workflow that generates relative
positions and orientations of two molecules on a regular trans-rot-rot grid,
and made the program, including ready-made scripts, available as the Python
package molgri.

The quality of a grid-based approach depends heavily on the selection of grid
points. To describe a binding process, we used a mesh of translation and
rotation grids. For translations, we employed a linear discretisation of the
radial distance between the two molecules. For rotations, we implemented six
different algorithms and compared them using a quantitative uniformity
measure. Algorithms based on polyhedra (Cube 3D and Icosahedron) yielded by
far the most uniform rotational grids. To us this was surprising, because
given that the quaternions most naturally describe rotations, we expected
quaternion-based algorithms to outperform the other algorithms.

By contrast, algorithms based on Euler angles missed some of the known minima
of the water dimer, even with high numbers of grid points, and should only be
used with caution for molecular systems. In terms of computational costs,
polyhedra-based algorithms showed the poorest scaling with increasing number
of grid points. However, since the run time to generate large grids with up to
1000 points is below 100 s on a standard workstation, the computational costs
for grid generation is negligible compared to the calculation of grid point
energies. Moreover, rotational grids are generated once and can be re-used for
any molecular system.

Grid points in state space are converted into corresponding molecular
structures and stored in a MD-trajectory file format, so that trans-rot-rot
grids can easily be interfaced with standard MD programs. The calculation of
the grid point energies can the thus be accomplished by already existing,
highly optimised programs. The entire method was first applied to two
approaching water molecules, showing that the position and depth of potential
energy minima found by molgri in combination with polyhedron- or quaternion-
based rotational grids matched those found by a stochastic sampling approach.
Moreover, we studied the radial and orientational dependency of energy
contributions, reproducing known oriented hydrogen bonding patterns for this
system and the dependencies of Lennard-Jones and Coulomb energies on
interaction distance. While this is not a new result for this simple system,
it demonstrates how the untangling of energy contributions on a grid could be
performed for more complex systems.

A drawback of our current approach is that it relies on the rigid body
approximation, which may be valid for small molecules but is insufficient for
induced fit interactions that are commonly encountered in macromolecular
systems. This problem can be tackled by replacing point energy calculations
with short energy minimisations. We demonstrated this on the system of a
coiled-coil peptide binding with a Cl- anion. The GROMACS energy minimzation
routine could handle starting structures with strongly overlapping Lennard-
Jones spheres, and in most cases the Cl- position could be accommodated by the
rotation of a side chain or a local adjustment of the backbone. This example
also showed that, at large distances, the long-range electrostatic
interactions are well-represented by the single-point energy calculations and
do not need an energy minimization. One can further improve the accuracy of
the grid point energies by running a short constraint simulation and by
including solvent molecules in the energy calculation. One could further
improve the accuracy of the grid point energies by running a short constraint
simulation and by including solvent molecules in the energy calculation.

We view the molgri package as a tool to represent the binding process on a
grid with approximate energies, to guide the choice of initial structures for
unconstrained sampling, or to guide the placement of umbrella potentials for
the construction of free-energy surfaces. With these use cases in mind, the
package can be extended in the following directions: ($i$) providing other
non-linear grids, e.g for cylindrical coordinates, ($ii$) improving the
interface with MD programs, such that molgri seamlessly calls a specific MD
program, ($iii$) extending the interface to quantum chemistry packages, such
that grid-point energies can be obtained from electronic-structure
calculations, ($iv$) extending the interface to flexible docking algorithms to
improve the handling of induced-fit binding. We are also interested in
extending the package to another use case: a rate-model of the binding process
constructed on a trans-rot-rot grid using the square-root approximation of the
Fokker-Planck equation [23, 24, 25, 26, 27]. For the square-root
approximation, the adjacency matrix of the grid, the high-dimensional volumes
of the grid cells, and the area of the interface between two neighboring grid
cells are required. Thus, full control over the grid is critical. On a
systematic trans-rot-rot grid, all mentioned parameters could be calculated
once and applied to various molecular systems.

Whether or not the grid-based approach could ever outperform sampling
approaches, assuming that ergodicity is fulfilled, is difficult to say. Two
aspects of the grid-based approach boost its computational efficiency: First,
because grid points are independent, calculation of the grid-point energies
can be completely parallelised. Second, because each grid point is evaluated
only once, a grid-based approach avoids re-examining the same minima over and
over again. Additionally, by their very nature, grid-based energy-calculations
should interface easily with adaptive resolution techniques [63]. On the other
hand, because of the high-dimensionality of conformational space, accurate
calculations of the grid point energies will have to rely on short stochastic
sampling of the orthogonal degrees of freedom in the foreseeable future like
we demonstrated in our minimisation routine. If the exploration of the
conformational space in orthogonal degrees of freedom is slow compared to the
degrees of freedom resolved by the grid, the computational efficiency is
dictated by the sampling.

However, the ergodicity assumption is a constant worry in MD simulations. In
contrast to sampling approaches, grid-based models of molecular interactions
do not need to assume ergodicity. Even with approximate estimates for the grid
point (free) energies, grid-based state exploration is a powerful complement
to sampling.

##  VII Author statements

###  VII.1 Competing interests statement

The authors declare there are no competing interests.

###  VII.2 Author contribution statement

H.Z.: Investigation; Methodology; Formal analysis; Software; Visualization;
Writing – original draft; Writing – review & editing. F.H.: Investigation;
Methodology; Formal analysis; Visualization; Writing – review & editing.
B.G.K.: Conceptualization; Methodology; Funding acquisition; Project
administration; Supervision; Writing – review & editing.

###  VII.3 Funding statement

This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through
grant SFB 1349 Fluorine-Specific Interactions - Project ID 387284271, Project
A05; SFB 1449 Dynamic Hydrogels at Biointerfaces - Project ID 431232613 -
Project C02; RTG 2473 Bioactive Peptides - Project ID 392923329 - Project C01.
We acknowledge access to high-performance computers via the Zentraleinrichtung
für Datenverarbeitung (ZEDAT) of Freie Universität Berlin.

###  VII.4 Data availability statement

Input molecular structure files used in this contribution and non-commercial
software used to manipulate them are included in our software package molgri.
We refer to Section III.3 for instructions on package acquisition and use.

##  VIII Supporting information

Detailed pseudo-algorithms and code dependencies for the software
implementation described in Section III are provided as supporting
information. Documentation and source code are hosted on GitHub
(https://github.com/bkellerlab/molecularRotationalGrids).

## References

  * Buckingham [1991] Buckingham, A. The hydrogen bond, and the structure and properties of H20 and (H20) 2. _Journal of molecular structure_ 1991, _250_ , 111–118.
  * Rinne et al. [2014] Rinne, K. F.; Gekle, S.; Netz, R. R. Ion-specific solvation water dynamics: single water versus collective water effects. _The Journal of Physical Chemistry A_ 2014, _118_ , 11667–11677.
  * Reinmuth et al. [2019] Reinmuth, M.; Pramanik, S.; Douglas, J. T.; Day, V. W.; Bowman-James, K. Structural impact of chelation on phytate, a highly phosphorylated biomolecule. _European Journal of Inorganic Chemistry_ 2019, _2019_ , 1870–1874.
  * Plattner and Noé [2015] Plattner, N.; Noé, F. Protein conformational plasticity and complex ligand-binding kinetics explored by atomistic simulations and Markov models. _Nature communications_ 2015, _6_ , 1–10.
  * Casasnovas et al. [2017] Casasnovas, R.; Limongelli, V.; Tiwary, P.; Carloni, P.; Parrinello, M. Unbinding kinetics of a p38 MAP kinase type II inhibitor from metadynamics simulations. _Journal of the American Chemical Society_ 2017, _139_ , 4780–4788.
  * Kahler et al. [2020] Kahler, U.; Kamenik, A. S.; Waibl, F.; Kraml, J.; Liedl, K. R. Protein-protein binding as a two-step mechanism: Preselection of encounter poses during the binding of BPTI and trypsin. _Biophysical journal_ 2020, _119_ , 652–666.
  * Wenz et al. [2021] Wenz, M. T.; Bertazzon, M.; Sticht, J.; Aleksic, S.; Gjorgjevikj, D.; Freund, C.; Keller, B. G. Target recognition in tandem WW domains: complex structures for parallel and antiparallel ligand orientation in h-FBP21 tandem WW. _Journal of Chemical Information and Modeling_ 2021,
  * Morriss-Andrews and Shea [2015] Morriss-Andrews, A.; Shea, J.-E. Computational studies of protein aggregation: methods and applications. _Annual review of physical chemistry_ 2015, _66_ , 643–666.
  * Morozova and Muthukumar [2018] Morozova, S.; Muthukumar, M. Electrostatic effects in collagen fibril formation. _The Journal of chemical physics_ 2018, _149_ , 163333\.
  * Hellmund et al. [2021] Hellmund, K. S.; Von Lospichl, B.; Böttcher, C.; Ludwig, K.; Keiderling, U.; Noirez, L.; Weiß, A.; Mikolajczak, D. J.; Gradzielski, M.; Koksch, B. Functionalized peptide hydrogels as tunable extracellular matrix mimics for biological applications. _Peptide Science_ 2021, _113_ , e24201.
  * Bolhuis et al. [2002] Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition Path Sampling: Throwing Ropes. _Annu. Rev. Phys. Chem_ 2002, _53_ , 291–318.
  * Barducci et al. [2011] Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2011, _1_ , 826–843.
  * Zuckerman and Chong [2017] Zuckerman, D. M.; Chong, L. T. Weighted ensemble simulation: review of methodology, applications, and software. _Annual review of biophysics_ 2017, _46_ , 43.
  * Bowman et al. [2010] Bowman, G. R.; Ensign, D. L.; Pande, V. S. Enhanced modeling via network theory: Adaptive sampling of Markov state models. _Journal of chemical theory and computation_ 2010, _6_ , 787–794.
  * Bruce et al. [2018] Bruce, N. J.; Ganotra, G. K.; Kokh, D. B.; Sadiq, S. K.; Wade, R. C. New approaches for computing ligand–receptor binding kinetics. _Current opinion in structural biology_ 2018, _49_ , 1–10.
  * Limongelli [2020] Limongelli, V. Ligand binding free energy and kinetics calculation in 2020. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2020, _10_ , e1455.
  * Plattner et al. [2017] Plattner, N.; Doerr, S.; De Fabritiis, G.; Noé, F. Complete protein–protein association kinetics in atomic detail revealed by molecular dynamics simulations and Markov modelling. _Nature chemistry_ 2017, _9_ , 1005–1011.
  * Capelli et al. [2019] Capelli, R.; Carloni, P.; Parrinello, M. Exhaustive search of ligand binding pathways via volume-based metadynamics. _The journal of physical chemistry letters_ 2019, _10_ , 3495–3499.
  * Lemke et al. [2018] Lemke, T.; Peter, C.; Kukharenko, O. Efficient Sampling and Characterization of Free Energy Landscapes of Ion–Peptide Systems. _Journal of chemical theory and computation_ 2018, _14_ , 5476–5488.
  * Batra et al. [2013] Batra, J.; Szabó, A.; Caulfield, T. R.; Soares, A. S.; Sahin-Tóth, M.; Radisky, E. S. Long-range electrostatic complementarity governs substrate recognition by human chymotrypsin C, a key regulator of digestive enzyme activation. _Journal of Biological Chemistry_ 2013, _288_ , 9848–9859.
  * Waldner et al. [2018] Waldner, B. J.; Kraml, J.; Kahler, U.; Spinn, A.; Schauperl, M.; Podewitz, M.; Fuchs, J. E.; Cruciani, G.; Liedl, K. R. Electrostatic recognition in substrate binding to serine proteases. _Journal of Molecular Recognition_ 2018, _31_ , e2727.
  * Chodera [2016] Chodera, J. D. A simple method for automated equilibration detection in molecular simulations. _Journal of chemical theory and computation_ 2016, _12_ , 1799–1805.
  * Bicout and Szabo [1998] Bicout, D.; Szabo, A. Electron transfer reaction dynamics in non-Debye solvents. _The Journal of chemical physics_ 1998, _109_ , 2325–2338.
  * Donati et al. [2018] Donati, L.; Heida, M.; Keller, B. G.; Weber, M. Estimation of the infinitesimal generator by square-root approximation. _Journal of Physics: Condensed Matter_ 2018, _30_ , 425201.
  * Heida [2018] Heida, M. Convergences of the squareroot approximation scheme to the Fokker–Planck operator. _Mathematical Models and Methods in Applied Sciences_ 2018, _28_ , 2599–2635.
  * Donati et al. [2021] Donati, L.; Weber, M.; Keller, B. G. Markov models from the square root approximation of the Fokker–Planck equation: calculating the grid-dependent flux. _Journal of Physics: Condensed Matter_ 2021, _33_ , 115902\.
  * Donati et al. [2022] Donati, L.; Weber, M.; Keller, B. G. A review of Girsanov Reweighting and of Square Root Approximation for building molecular Markov State Models. _arXiv preprint arXiv:2209.10544_ 2022,
  * Euler [1776] Euler, L. Nova methodus motum corporum rigidorum degerminandi. _Novi commentarii academiae scientiarum Petropolitanae_ 1776, 208–238.
  * Kirk [2012] Kirk, D. _Graphics Gems III (IBM Version)_ ; Elsevier, 2012.
  * Yershova and LaValle [2004] Yershova, A.; LaValle, S. M. Deterministic sampling methods for spheres and SO (3). IEEE International Conference on Robotics and Automation. 2004; pp 3974–3980.
  * Yershova et al. [2010] Yershova, A.; Jain, S.; Lavalle, S. M.; Mitchell, J. C. Generating uniform incremental grids on SO (3) using the Hopf fibration. _The International journal of robotics research_ 2010, _29_ , 801–812.
  * Karney [2007] Karney, C. F. Quaternions in molecular modeling. _Journal of Molecular Graphics and Modelling_ 2007, _25_ , 595–604.
  * Diebel [2006] Diebel, J. Representing attitude: Euler angles, unit quaternions, and rotation vectors. _Matrix_ 2006, _58_ , 1–35.
  * Purser and Rančić [2011] Purser, R. J.; Rančić, M. _A standardized procedure for the derivation of smooth and partially overset grids on the sphere, associated with polyhedra that admit regular griddings of their surfaces; Mathematical principles of classification and construction; Part I_ ; 2011.
  * Sadourny [1972] Sadourny, R. Conservative finite-difference approximations of the primitive equations on quasi-uniform spherical grids. _Monthly Weather Review_ 1972, _100_ , 136–144.
  * Lindemann et al. [2004] Lindemann, S. R.; Yershova, A.; LaValle, S. M. _Algorithmic Foundations of Robotics VI_ ; Springer, 2004; pp 313–328.
  * Kneller [1991] Kneller, G. R. Superposition of molecular structures using quaternions. _Molecular Simulation_ 1991, _7_ , 113–119.
  * Fincham [1992] Fincham, D. Leapfrog rotational algorithms. _Molecular Simulation_ 1992, _8_ , 165–178.
  * Kol et al. [1997] Kol, A.; Laird, B. B.; Leimkuhler, B. J. A symplectic method for rigid-body molecular simulation. _The Journal of chemical physics_ 1997, _107_ , 2580–2588.
  * Miller Iii et al. [2002] Miller Iii, T.; Eleftheriou, M.; Pattnaik, P.; Ndirango, A.; Newns, D.; Martyna, G. Symplectic quaternion scheme for biophysical molecular dynamics. _The Journal of chemical physics_ 2002, _116_ , 8649–8659.
  * Nielsen et al. [2010] Nielsen, S. O.; Bulo, R. E.; Moore, P. B.; Ensing, B. Recent progress in adaptive multiscale molecular dynamics simulations of soft matter. _Physical Chemistry Chemical Physics_ 2010, _12_ , 12401–12414.
  * Stumpe and Grubmüller [2007] Stumpe, M. C.; Grubmüller, H. Aqueous urea solutions: structure, energetics, and urea aggregation. _The Journal of Physical Chemistry B_ 2007, _111_ , 6220–6228.
  * Heinz and Grubmüller [2019] Heinz, L. P.; Grubmüller, H. Computing spatially resolved rotational hydration entropies from atomistic simulations. _Journal of Chemical Theory and Computation_ 2019, _16_ , 108–118.
  * Lynden-Bell and Stone [1989] Lynden-Bell, R.; Stone, A. Reorientational correlation functions, quaternions and Wigner rotation matrices. _Molecular Simulation_ 1989, _3_ , 271–281.
  * Wikipedia contributors [2022] Wikipedia contributors, Rotation matrix - Wikipedia, The Free Encyclopedia. 2022; https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions, [accessed 24-October-2022].
  * Hamilton [1840] Hamilton, W. R. On a new species of imaginary quantities, connected with the theory of quaternions. _Proceedings of the Royal Irish Academy (1836-1869)_ 1840, _2_ , 424–434.
  * Hamilton [1850] Hamilton, W. R. On quaternions and the rotation of a solid body. Proceedings of the Royal Irish Academy. 1850; pp 38–56.
  * Evlero and Petropolitanae [1770] Evlero, A. L.; Petropolitanae, N. Problema algebraicvm ob affectiones prorsvs singvlares memorabile. _Novi Commentarii Academiae imperialis scientiarum Petropolitanae_ 1770,
  * Rodrigues [1840] Rodrigues, O. Des lois géométriques qui régissent les déplacements d’un système solide dans l’espace, et de la variation des coordonnées provenant de ces déplacements considérés indépendamment des causes qui peuvent les produire. _J. Math. Pures Appl_ 1840, _5_ , 5.
  * Lynch and Park [2017] Lynch, K. M.; Park, F. C. _Modern robotics_ ; Cambridge University Press, 2017\.
  * Wikipedia contributors [2022] Wikipedia contributors, Spherical cap — Wikipedia, The Free Encyclopedia. 2022; {https://en.wikipedia.org/wiki/Spherical_cap}, [accessed 24-October-2022].
  * mol [2022] pypi package molgri. 2022; https://pypi.org/project/molgri/, [accessed 29-October-2022].
  * mol [2022] GitHub repository molgri. 2022; https://github.com/bkellerlab/molecularRotationalGrids, [accessed 29-October-2022].
  * Woolfson [2005] Woolfson, D. N. The design of coiled-coil structures and assemblies. _Advances in protein chemistry_ 2005, _70_ , 79–112.
  * Van Der Spoel et al. [2005] Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. GROMACS: fast, flexible, and free. _Journal of computational chemistry_ 2005, _26_ , 1701–1718.
  * Abraham et al. [2015] Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. _SoftwareX_ 2015, _1_ , 19–25.
  * Lindahl et al. [2001] Lindahl, E.; Hess, B.; Van Der Spoel, D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. _Molecular modeling annual_ 2001, _7_ , 306–317.
  * Berendsen et al. [1995] Berendsen, H. J.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. _Computer physics communications_ 1995, _91_ , 43–56.
  * Humphrey et al. [1996] Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. _Journal of molecular graphics_ 1996, _14_ , 33–38.
  * Harrington and Engelen [2004] Harrington, B.; Engelen, J. Inkscape. 2004; Software available at http://www. inkscape. org.
  * Jorgensen et al. [1983] Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. _The Journal of chemical physics_ 1983, _79_ , 926–935.
  * Lindorff-Larsen et al. [2010] Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. _Proteins: Structure, Function, and Bioinformatics_ 2010, _78_ , 1950–1958.
  * Cortes-Huerto et al. [2021] Cortes-Huerto, R.; Praprotnik, M.; Kremer, K.; Delle Site, L. From adaptive resolution to molecular dynamics of open systems. _The European Physical Journal B_ 2021, _94_ , 1–22.