††thanks: present address: Schrödinger Inc., 1540 Broadway 24th Floor, New
York, NY 10036, United States

# Systematic Analysis of Biomolecular Conformational Ensembles with PENSA

Martin Vögele  Department of Computer Science, Stanford University School of
Engineering  Department of Molecular and Cellular Physiology, Stanford
University School of Medicine  Department of Structural Biology, Stanford
University School of Medicine  Neil J. Thomson  Division of Computational
Biology, School of Life Sciences, University of Dundee  Sang T. Truong
Department of Computer Science, Stanford University School of Engineering
Jasper McAvity  Department of Computer Science, Stanford University School of
Engineering  Ulrich Zachariae  Division of Computational Biology, School of
Life Sciences, University of Dundee  Ron O. Dror  Department of Computer
Science, Stanford University School of Engineering  Department of Molecular
and Cellular Physiology, Stanford University School of Medicine  Department of
Structural Biology, Stanford University School of Medicine  Institute for
Computational and Mathematical Engineering, Stanford University

###### Abstract

Molecular simulations enable the study of biomolecules and their dynamics on
an atomistic scale. A common task is to compare several simulation conditions
- like mutations or different ligands - to find significant differences and
interrelations between them. However, the large amount of data produced for
ever larger and more complex systems often renders it difficult to identify
the structural features that are relevant for a particular phenomenon. We
present a flexible software package named PENSA that enables a comprehensive
and thorough investigation into biomolecular conformational ensembles. It
provides a wide variety of featurizations and feature transformations that
allow for a complete representation of biomolecules like proteins and nucleic
acids, including water and ion cavities within the biomolecular structure,
thus avoiding bias that would come with manual selection of features. PENSA
implements various methods to systematically compare the distributions of
these features across ensembles to find the significant differences between
them and identify regions of interest. It also includes a novel approach to
quantify the state-specific information between two regions of a biomolecule
which allows, e.g., the tracing of information flow to identify signaling
pathways. PENSA also comes with convenient tools for loading data and
visualizing results in ways that make them quick to process and easy to
interpret. PENSA is an open-source Python library maintained at
https://github.com/drorlab/pensa along with an example workflow and a
tutorial. Here we demonstrate its usefulness in real-world examples by showing
how it helps to determine molecular mechanisms efficiently.

††preprint: AIP/123-QED

##  I Introduction

Molecules exist not as static structures, but in a range of conformations that
fluctuate about energetic equilibria and can be described as a thermodynamic
ensemble. In recent years, molecular dynamics (MD) simulations have become one
of the standard methods in molecular biology,Hollingsworth and Dror (2018)
providing detailed insights into a molecule’s conformational ensemble,
complementing static experimental structures that represent only the most
probable conformation.Knoverek  _et al._ (2018); Suomivuori  _et al._ (2020) A
common and important problem in molecular biophysics is to probe the effect of
two or more similar conditions on a macromolecule and identify the resulting
conformational differences and how they relate to each other. Typical analyses
include investigating the effects of: small-molecule ligand binding Provasi
_et al._ (2011); McCorvy  _et al._ (2018), protein mutations, Cordero-Morales
_et al._ (2007) or protonation changesLiu  _et al._ (2015); Thomson  _et al._
(2021), usually by simulating the macromolecule with and without a particular
ligand, mutation, or protonation state. Visual inspection is commonly
performed to identify the resulting differences and to derive causal
interrelations.Humphrey  _et al._ (1996); DeLano (2002) However, in systems
where the subtleties of a small shift in the populations of minority
conformational states or the motion of single atoms might have crucial
functional outcomes, analysis is complicated by the large dimensionality of
most systems. Hollingsworth and Dror (2018) Even when the dimensionality is
distilled to a subset of features, such as torsion angles, simulations contain
billions of frames, and functionally relevant differences can extend to
regions far from a small-molecule binding site or mutation site.Provasi  _et
al._ (2011); Dror  _et al._ (2011); Zhou  _et al._ (2019); Dror  _et al._
(2013); Bowman  _et al._ (2015); Hollingsworth and Dror (2018) The necessary
analysis is the bottleneck of many molecular biology projects as it can take
weeks of dedicated work if performed by eye and by one-off scripts, and a
focus on preconceived candidate mechanisms can lead to missing unexpected
effects.

In light of these hurdles, the strong interest in ensemble analyses over the
past two decades has led to development of ensemble databases with inbuilt
analysis tools Zivanovic  _et al._ (2020) and the availability of more
powerful, systematic and quantitative approaches, including: single-score
similarity measures between two ensembles,Brüschweiler (2003); Lindorff-Larsen
and Ferkinghoff-Borg (2009), implemented in libraries such as EncoreTiberti
_et al._ (2015); generative deep learning methods;Noé  _et al._ (2019) Markov
model analyses; Husic and Pande (2018); Nüske  _et al._ (2017) and neural-
network based analyses.Fraccalvieri  _et al._ (2011); Ward  _et al._ (2021)
For example, DiffNetsWard  _et al._ (2021) successfully identified mutation
sites that affect the signalling profile of the oxytocin receptor.Malik  _et
al._ (2021) However, the available methods are either computationally costly,
difficult to apply, have to be fine-tuned for a particular system, or are not
easily interpretable. This shortage of standardized yet flexible analysis
tools, specifically for comparison of multiple conditions, poses an obstacle
in many research avenues.

We present the modular software library PENSAVögele  _et al._ (2021) (short
for Python ENSemble Analysis) that enables the flexible implementation of
systematic and quantitative yet easily interpretable workflows for exploratory
analysis of biomolecular conformational ensembles (Fig. 1). It contains user-
friendly tools to preprocess simulation data, to apply various analysis
methods across simulation conditions, and to visualize the results. PENSA is
an open-source Python library maintained at https://github.com/drorlab/pensa
along with an example workflow and a tutorial.

PENSA reduces the large dimensionality of molecular systems to features that
allow for its complete representation (Fig. 2). A typical PENSA workflow first
determines the same features for all ensembles, currently including: amino
acid torsions, nucleic acid torsions, arbitrary interatomic distances, and a
novel featurization method for water cavities and ion binding sites. The
library provides the necessary feature readers and preprocessing tools while
allowing for the integration of custom implementations. The ensemble
definition is flexible, e.g., the data representing two ensembles could be
derived from two independent simulations or from two parts of a single
simulation that encompasses two ensembles. Taking into account all features,
PENSA attenuates the bias that would come with manual pre-selection. To
further reduce the complexity of a system, classical methods for
dimensionality-reduction and clustering can be applied across the joint
ensemble combining all conditions, and multiple primary features can be
combined to one via multivariate discretization. These features form the basis
for the subsequent quantitative analysis that provides a comprehensive insight
into the ensembles (Fig. 1).

Figure 1: PENSA allows a variety of analysis workflows. Green boxes show
primary features (section II.1) and red boxes show possible feature
transformations (section II.2). The blue box represents comparison using
either the Jensen–Shannon distance (JSD) or the Kolmogorov Smirnov Statistic
(KSS), and the yellow box represents the quantification of mutual information
using State-Specific Information (SSI) or State-Specific Co-information
(CoSSI), which are discussed in sections II.3 and II.4, respectively. PENSA’s
modular nature allows for flexibility in the path from input data to analysis
metric.

PENSA’s analysis methods focus on quantitatively exploring multiple
conformational ensembles and discovering interrelations within and between
them. To locate the most relevant differences between two ensembles, PENSA
includes a direct comparison of every feature’s distribution between ensembles
via Jenson-Shannon distance (JSD) and the Kolmogorov-Smirnov statistic (KSS).
Furthermore, a mutual information analysis based on state-specific information
Thomson  _et al._ (2021) (SSI), is included to provide a measure of the
information that features signal about the ensembles’ conditions or the
transitions between them. An extension of SSI to three variables - CoSSI,
enables the tracing of information flow between two regions along the
intermediary features. While mostly designed with the comparison of two
simulation conditions in mind, PENSA’s analysis methods can be further
expanded to more than two ensembles. For an accessible interpretation of the
analysis, various visualization options exist that conveniently transfer
results to plots, heat maps, or project results onto three dimensional
reference structures. With these methods, PENSA makes it easy to perform
quantitative and systematic exploratory analyses across multiple conditions
for a variety of biomolecular systems.

Figure 2: Examples for biomolecular features implemented in PENSA. _A_ : Amino
acid torsions. $\Phi$ and $\Psi$ represent amino acid backbone torsion angles,
while $\chi_{i}$ represent amino acid side chain torsion angles, as depicted
on Tryptophan. _B_ : Representative angles for DNA and RNA. Six backbone-
torsions sufficiently represent the torsional modes of DNA and RNA. _C_ :
Distances between backbone C-$\alpha$’s of residue pairs in two G-protein
coupled receptor ensembles. These distances represent divergences between
ensembles in example study seen in section III.3. _D_ : Description of top 10
most probable water sites featurized from water density surrounding the
protein accompanied with close up of featurization. Structures are visualized
using PyMolSchrödinger, LLC (2015) and density features using Biotite Kunzmann
and Hamacher (2018).

Here, we discuss the functionality included in PENSA and demonstrate its
usefulness on three real-world applications: We show how to describe the
influence of local frustration on loop opening during the catalytic cycle of
an oxidoreductase Stelzl  _et al._ (2020), how to quantify the influence of
force-field parameter changes on simulations of nucleic acids Cruz-León  _et
al._ (2021), and how to identify the effect of ionization on receptor
proteins. The examples mostly reproduce existing results, confirming the
reliability of our approach, but also show additional new discoveries made
possible by systematic comparison of molecular simulations.

##  II PENSA Functionality

###  II.1 Featurization

In general, the 3D structure of a molecule is represented with Cartesian
coordinates, but in molecular dynamics, this representation is mostly
redundant due to bonded and non-bonded atomic constraints. Other
representations summarise the motion of multiple atomic coordinates in a
succinct and instructive manner. In the following, we use the term feature for
every such numerical attribute calculated from the coordinates of a structure.
PENSA represents each ensemble (e.g., simulation trajectory) by the same set
of $F$ features $\\{x_{1},...,x_{F}\\}$. Here we discuss the types of features
we consider most suited for a systematic analysis of molecular structure,
currently including torsion angles, interatomic distances, water cavities and
ion binding sites (for an overview, see Fig. 2). They are implemented in PENSA
using MDAnalysis Michaud-Agrawal  _et al._ (2011); Gowers  _et al._ (2016) and
PyEMMA Scherer  _et al._ (2015) and can be replaced or extended by other
descriptions. PENSA is modular and new feature readers can be added easily.

_Representing Molecules by Torsion Angles:_ Because the energy minima of
molecular structures are largely determined by stereochemistry, the 3D
conformation of most biological macromolecules can be approximately described
using a set of torsion angles. For proteins, we use torsions around chemical
bonds to describe the effective degrees of freedom in the protein backbone and
the side chains (Figure 2). The three rotatable bonds in the protein backbone
are characterized by angles $\phi$, $\psi$, and $\omega$. As $\omega$ is
almost always in trans configuration, we only use $\phi$ and $\psi$. In
addition, the amino acid side chains have up to five rotatable bonds. This
results in a number of variables that grows linearly with the number of
residues, describing the full protein structure while avoiding the redundancy
that comes with using Cartesian coodinates (Fig. 2A). The analogous
description of nucleic acids via torsion angles requires the definition of
pseudo-torsionsKeating  _et al._ (2011). The 3D structure of DNA or RNA can be
described via the six main chain torsion angles ($\alpha$, $\beta$, $\gamma$,
$\delta$, $\epsilon$, $\zeta$) around the covalent bonds and the angle $\chi$
about the glycosidic bond and the sugar pucker (Fig. 2B). An alternative
description — also implemented in PENSA — uses two pseudo-torsions, $\eta$ and
$\theta$, that are defined around imaginary lines connecting more distant
atoms.Keating  _et al._ (2011)

_Representing Molecules by Distance:_ The set of all distances between all
atoms provides a complete description of the structure of a molecule,
independent of the coordinate system. However, taking into account all $N$
atoms is extremely inefficient as the number of distances required grows
$\propto N^{2}$ while the number of degrees of freedom grows $\propto N$.
Reducing the set of distances to a few relevant ones can make this approach
much more efficient. For example, the overall structure of a protein is
usually described by the distances between all C$\alpha$ atoms. A further
reduction in redundancy can be achieved by using system-specific knowledge,
obtaining a subset that is smaller but still representative for the overall
dynamics of the system. Common examples for this are hydrogen bond lengths to
describe bonding patterns or base-pair distances in DNA. A combined analysis
of distances and torsions can be beneficial as some effects can be easier to
spot in distances and others in torsions.

_Representing Water Cavities:_ We have implemented a new method to represent
the presence and orientation of internal water molecules, which are crucial to
the structure and function of many biomolecules. Typically, internal water
molecules are defined by the protein sites they bind to, or by an atomic
density that averages the motion of the waters.Venkatakrishnan  _et al._
(2019); Yuan  _et al._ (2014); Pardo  _et al._ (2007); Levy and Onuchic (2004)
PENSA adds to this with a dynamic, orientation-based representation of water
molecules in cavities within a biomolecule. Because water cavities are often
accessible to water molecules in a freely diffusing bulk solvent, water
molecules may interchange in a water cavity while maintaining the function of
the cavity. Therefore, featurization of a certain water _site_ within the
biomolecule is more important than the individual water molecules that occupy
the cavity. The polar nature of water molecules enables a water cavity to
function as a polarizable interface, mediating hydrogen bond networks between
amino acids that would otherwise be impossible to form. Any cavity that can
accommodate a water molecule can additionally be unoccupied, whereby the
occupation of the cavity can act as a further feature. To account for both of
these effects, PENSA featurizes a water cavity via its occupation and, when
occupied, the orientation of the water molecule. First, we locate the cavities
by finding all local maxima within a 3D density grid. The grid is obtained by
aligning the assembled trajectories such that the density refers to the
combined ensemble and can be used to compare identical cavities across the
ensembles. Each local maximum marks the center of a water cavity, defined as a
sphere of radius 3.5 Å, based on hydrogen bonding interaction
distances.Kajander  _et al._ (2000) For each simulation frame, the cavity
feature value is defined by the orientation of the water molecule (if
occupied), or a value representing its unoccupied state. The water molecule’s
orientation is represented by the angular components of its dipole moment in
spherical coordinates. This representation is appropriate when rotation,
translation, and periodicity effects of the simulation system during the
simulation are excluded from the water pocket featurization during
preprocessing.

_Representing Ion Cavities:_ PENSA provides a similar approach to featurize
ion cavities, which play important roles in the structure and function of many
biomolecular systems — from ion channels to catalytic enzymes Zarzycka  _et
al._ (2019); Andreini  _et al._ (2008). Simulation trajectories must be
preprocessed in the same manner as for water cavities, ensuring that the ion
density grid is not affected by any system rotation, translation or switches
between periodic boundaries. Ion cavities are then represented as a binary
feature describing the presence or absence of an ion in the cavity, or as a
discrete feature describing the ion indices that occupy the cavity, for
instance in order to keep track of functionally relevant ion binding and
unbinding events. This representation allows for an investigation of the
movement of specific ions in and out of an array of ion binding sites, and has
already been employed to identify co-operative knock-on permeation of ions
bound to different sites of cation channels.Ives  _et al._ (2022)

###  II.2 Feature Transformations

_Multivariate Discretization:_ PENSA offers an automatic discretization of one
or multiple combined features into distinct states. For example, amino acid
torsions typically oscillate about local minima in populations known as
rotamers, often sampling different rotamers in an MD simulation.Dunbrack
(2002); Dunbrack and Karplus (1993); Scouras and Daggett (2011) In the
Dunbrack library,Dunbrack and Karplus (1993), each local rotamer can be
accurately represented by a Gaussian distribution around the minimum energy
conformation. Similarly, Gaussians can represent the oscillation of a water
molecule’s polarization (orientation of its dipole) within a water cavity,
while cavity occupancy changes can be conveniently considered as two discrete
populations with zero oscillation. The discrete states of these distributions
are generated by applying a multi-modal Gaussian fit. The fit parameters are
obtained with a non-linear least squares fit of up to ten Gaussians to each
feature distribution using the SciPy library Virtanen  _et al._ (2020). For
our purposes, we found this method to be computationally more efficient than
alternative methods such as Gaussian Mixture Model. The state limits for each
distribution are then defined by the Gaussian intersects. Distributions which
have a cyclic periodicity may oscillate about a periodic boundary but show up
as two states on different ends of the distribution. To account for this
boundary effect, all periodic distributions can be linearly shifted. In
general, state limits can be defined for any kind of distribution, even non-
Gaussian, using suitable clustering algorithms. For these cases, PENSA allows
the manual input of arbitrary state boundaries. In instances where the
dynamics of an amino acid are described by a combined view of all side-chain
angles, e.g., the five side-chain torsions of arginine, an arbitrary number of
N dimensions can be combined into one joint feature. The states of an
N-dimensional feature are then the combination of all discretized microstates,
defining a grid in the combined feature space, as shown in Fig. 3.

Figure 3: Discrete macrostates of a water cavity generated by discretization
of individual pairs of spherical coordinates $(\phi,\theta)$ of the polar
vector of each water molecule within the water cavity. Red dashed lines depict
the state boundaries and include the periodic boundaries as the first and last
state boundaries.

_Dimensionality Reduction:_ Comparisons of multiple complex systems can be
made more tractable by representing these systems in terms of low-dimensional,
higher-order descriptions that summarize a large number of primary features in
a small number of quantities. The most popular algorithm for such
dimensionality reduction is Principal Component Analysis (PCA) F.R.S. (1901).
Time series, such as those from molecular dynamics simulations, can also be
projected via Time-Lagged Independent Component Analysis (TICA) Schultze and
Grubmüller (2021). Classically, dimensionality reduction has been applied on
single ensembles. But to detect patterns across ensembles and differences
between them, we have to define the same representation for all investigated
ensembles. PENSA users can perform dimensionality reduction on the combined
data of all ensembles included in the analysis and then compare them along the
resulting reduced dimensions.

_Clustering:_ The coordinate space defined by the dimensionality reduction
methods discussed above can be clustered into discrete states, again using
data from all ensembles. Clustering the structures from all ensembles in the
resulting lower-dimensional space provides discrete states. PENSA implements
k-means clustering Macqueen (1967) and regular-space clustering Prinz  _et
al._ (2011), two popular algorithms for this task. Users can calculate
populations of the resulting discrete states in each ensemble and compare
them.

###  II.3 Feature-by-Feature Comparison

The local extent of deviations between two ensembles can be quantified by
comparing each feature’s probability distribution in one ensemble to its
distribution in the other ensemble. For a feature $x_{f}$ in a simulation
trajectory of length $T$, the $T$ samples of $x_{f}$ give an empirical
estimation of the distribution $p(x_{f})$, which describes the behaviour of
$x_{f}$ in that ensemble. We want to compare the distribution $p_{i}(x_{f})$
in ensemble $i$ to the corresponding distribution $p_{j}(x_{f})$ in ensemble
$j$. These distributions may have complex functional form, hence comparing
their summary statistics (mean and standard deviation) might not be sufficient
as they cannot capture more subtle differences in the distributions, e.g., the
split-up of one state (unimodal distribution) into two states (bimodal).
Instead, PENSA provides comparison measures that are designed to capture
differences in probability distributions, namely the Jensen-Shannon Distance
(JSD) and the Kolmogorov-Smirnov Statistic (KSS) divergence.

_Jensen-Shannon Distance:_ Two distributions can be compared using the Jensen-
Shannon distance $D_{JS}$, a symmetrized and numerically more stable version
of the Kullback-Leibler divergence $D_{KL}$. For two distributions over a
feature $x_{f}$ from ensemble $i$ and ensemble $j$ , $D_{JS}$ is defined below

$\displaystyle D_{\mathrm{JS}}[p_{i}\|p_{j}]$ $\displaystyle=$
$\displaystyle\frac{1}{4}\sqrt{D_{KL}[p_{j}\|\bar{p}]+D_{KL}[p_{j}\|\bar{p}]}$
(1) $\displaystyle\mathrm{with}\;\bar{p}$ $\displaystyle=$
$\displaystyle\frac{p_{i}+p_{j}}{2}$

with $D_{KL}$ the Kullback-Leibler divergence. For numerical reasons, we
always use its discrete version:

$\displaystyle D_{KL}[p_{i}(x_{f})||p_{j}(x_{f})]$ $\displaystyle=\int
p_{j}(x_{f})\log\frac{p_{i}(x_{f})}{p_{j}(x_{f})}dx$ (2)
$\displaystyle\approx\sum_{x_{f}\in\mathcal{X}_{f}}p_{j}(x_{f})\log\frac{p_{i}(x_{f})}{p_{j}(x_{f})}$

where ${\mathcal{X}}_{f}$ is the set of possible states. In the case of
continuous features, these states are bins along the feature coordinate,
obtained by evenly dividing the range of the joint distribution. Note that
$D_{KL}$ is not symmetric, i.e., $D_{\mathrm{KL}}[p_{i}\parallel p_{j})\neq
D_{\mathrm{KL}}[p_{j}\parallel p_{i}]$ but $D_{JS}$ is. The use of JSD as a
comparison metric has been discussed in more detail in previous workLindorff-
Larsen and Ferkinghoff-Borg (2009) where the comparison was performed on
entire ensembles instead of individual features. Most importantly, in contrast
to the unbounded — and in practice often divergent — KL divergence, the
Jensen-Shannon distance ranges from 0 to 1 where 0 is obtained for identical
distributions and 1 for a pair of completely different distributions.

_Kolmogorov-Smirnov Statistic:_ Alternatively, we can quantify the differences
between two distributions of continuous features without the need to define a
binning parameter by using the Kolmogorov–Smirnov statistic $D_{KS}$. It is
defined for a feature $x_{f}$ as

$D_{\mathrm{KS}}(p_{i}\|p_{j})=\sup_{x_{f}}\|F_{p_{i}}(x_{f})-F_{p_{j}}(x_{f})\|$
(3)

with $F_{p_{i}}(x_{f})$ and $F_{p_{j}}(x_{f})$ the empirical distribution
functions of $p_{i}(x_{f})$ and $p_{f}(x_{f})$, respectively. The empirical
distribution functions are directly obtained from the calculated features and
require no sorting of the data into arbitrary bins. The results of
$D_{\mathrm{JS}}$ and $D_{\mathrm{KS}}$ for the same comparison ideally are
very similar which can serve as an important sanity check.

_Overall Ensemble Similarity:_ The overall similarity of two ensembles over
all features in a metrics can be quantified by aggregating similarity metrics
of all features $x_{f}$. For example, an average Kolmogorov-Smirnov statistic
$\bar{D}_{KS}$ of two ensembles $i$ and $j$ can be computed as:

$\bar{D}_{KS}=\frac{1}{F}\sum_{f=1}^{F}D_{KS}[p_{i}(x_{f})||p_{j}(x_{f})]$ (4)

and $\bar{D}_{JS}$ analogously. Aggregation functions other than the average,
including the maximum and the minimum, are also implemented. Similar to other
metrics that quantify the similarity of two ensembles in a single
score,Brüschweiler (2003); Lindorff-Larsen and Ferkinghoff-Borg (2009);
Tiberti  _et al._ (2015) these aggregated metrics are particularly helpful
when we need to evaluate the output of a new method to a reference ensemble,
for example comparing a simulation or a generative machine learning model to a
ground truth from an experiment or a more accurate level of simulation.

###  II.4 Mutual Information Analysis

A mutual information analysis can be employed to measure how much the specific
value of one feature is coupled to the specific value of another.Shannon
(1948) Applying this approach to two conformational ensembles $i$ and $j$, one
can identify if the specific values of a feature are dependent on the system
ensemble – $i$ or $j$, and vice versa, how much the value of a feature reveals
about whether it stems from ensemble $i$ or $j$. The ensemble identifiers $i$
and $j$ are then equivalent to the values of an additional feature within one
(joint) ensemble. PENSA focuses on mutual information shared between a
feature’s conformational states (e.g., the multivariate states described in
section II.2) and the ensemble states $i$ and $j$. To quantify this, we
calculate the State-Specific Information (SSI, Fig. 4), a linear, discrete-
state adaptation of mutual information that has originally been developed for
amino acid torsions that act as molecular switchesThomson  _et al._ (2021) but
here is generalized to any feature with a distribution that can be represented
as discrete states. It should be noted that an arbitrary number of ensembles
could be incorporated into the SSI calculation, simultaneously measuring the
mutual information between all conditions and a feature, but it is currently
implemented for two. Similarly, SSI can operate on a single ensemble that is
partitioned into two sub-ensembles, e.g., along a state boundary.

Figure 4: Schematic illustration of the calculation process for the State-
Specific Information of feature $x_{f}$, $I_{SSI}(x_{f})$. The combined data
points for feature $x_{f}$ are discretized into states defined by the
intersects of a Gaussian fit to the data histogram; blue and orange represent
ensembles i and j, respectively. Each data point is transformed into the
corresponding discrete state it lies within ($S(x_{f})$) and combined with the
ensemble states - $i$ or $j$. Colors are added to the data points on step 4 to
highlight the unique substates that arise from the combination of feature
states and ensemble states.

_State-Specific Information (SSI):_ The SSI measure $I_{\mathrm{SSI}}(x_{f})$
quantifies the degree to which conformational state transitions of feature
$x_{f}$ signal information about the ensembles $i$ and $j$ or the transitions
between them. It is defined as

$I_{\mathrm{SSI}}(x_{f})=\sum\limits_{\begin{subarray}{c}s\in S(x_{f})\\\
e\in[i,j]\end{subarray}}p(s,e)\log\frac{p(s,e)}{p(s)p(e)}$ (5)

where $S(x_{f})$ is the transformation of the combined data points of $x_{f}$
across ensembles $i$ and $j$ into data points referring to a feature state-
identifier $s$, derived from the discretized probability distribution of
$x_{f}$, as seen in figure 4. The variable $e$ is a generic representation of
the ensemble ID, with possible values $i$ or $j$. $p(s,e)$ is the joint
distribution of a feature state of $x_{f}$ and ensemble state (obtained from
the probabilities of each element in the list shown in Fig. 4, step 4.), and
$p(s)$ and $p(e)$ are the marginal distributions. The SSI ranges from 0 bits
to 1 bit, where 0 bits represents no shared information and 1 bit represents
maximal shared information between the ensemble (transitions) and the
features.

_State-Specific Co-Information (CoSSI):_ To quantify the degree to which two
features interact with one another as they signal information about the
ensemble they are in ($i$ or $j$), State-Specific Co-Information (CoSSI) is
employed. This multivariate feature-feature-ensemble metric is the linear,
discrete-state adaptation of co-informationBell (2003) that uses Shannon’s
discrete entropy formulation Shannon (1948), calculated using

$\displaystyle I_{\mathrm{CoSSI}}(x_{1},x_{2})$ (6)
$\displaystyle=\sum\limits_{\begin{subarray}{c}s_{1}\in S(x_{1})\\\ s_{2}\in
S(x_{2})\end{subarray}}p(s_{1},s_{2})\log\frac{p(s_{1},s_{2})}{p(s_{1})p(s_{2})}$
$\displaystyle-\sum\limits_{\begin{subarray}{c}s_{1}\in S(x_{1})\\\ s_{2}\in
S(x_{2})\\\
e\in[i,j])\end{subarray}}p(s_{1},s_{2},e)\log\frac{p(s_{1},s_{2},e)p(e)}{p(s_{1},e)p(s_{2},e)}$

where the transformation $S(x_{f})$ is as previously defined. $I_{CoSSI}$ can
be positive or negative, indicating whether the switch between ensembles
increases ($I_{CoSSI}>0$), decreases ($I_{CoSSI}<0$), or does not affect
($I_{CoSSI}=0$) the communication between two features $x_{1}$ and $x_{2}$. In
the case of small-molecule ligand binding, for instance, positive $I_{CoSSI}$
between features can represent the turning-on of a signal channel by a ligand.

###  II.5 Visualization

PENSA includes convenient functions to visualize all stages of the analysis
workflow. Primary features as well as processed features (like projections
onto PCA eigenvectors) can easily be compared individually using histograms or
inferred densities, and combinations of two features using heatmaps, with all
functionality based on Matplotlib.Hunter (2007) The water and ion cavity
featurizer writes out a structure file with the average position of the
molecular ensemble, generated via MDAnalysis,Michaud-Agrawal  _et al._ (2011)
with additional atoms added via BiotiteKunzmann and Hamacher (2018) that mark
the cavity centers and store the magnitude of the probability maxima. Input
structures can be sorted along the values of primary or processed features
using MDAnalysisMichaud-Agrawal  _et al._ (2011) which is particularly useful
for PCA or tICA to see which component of a molecule’s motion is associated
with which eigenvector. Analysis metrics for comparison and mutual information
that are related to a single residue (e.g., the maximum JSD of all side-chain
torsions in an amino acid) can be stored in structure files using
MDAnalysisMichaud-Agrawal  _et al._ (2011) and we provide scripts for
PyMolSchrödinger, LLC (2015) and VMDHumphrey  _et al._ (1996) to visualize
them via the color or the width of the cartoon representation. Metrics related
to two features (e.g., distances) are visualized in square heatmaps, also
implemented via Matplotlib.Hunter (2007) These visualization options provide a
comprehensive overview of complex systems in one — or very few — figures.

##  III Example Applications

###  III.1 Understanding effects of a small chemical modification: loop
opening in an oxidoreductase

As a first example, we show how the systematic comparison of protein backbone
and side-chain torsions provides a comprehensive overview on the differences
in the conformational ensembles induced by a small chemical modification. We
consider the oxidation of two cysteine thiols to a disulfide bond in the
N-terminal domain of the key bacterial oxidoreductase DsbD (nDsbD). DsbD plays
an important role in electron transport across the inner cytoplasmic membrane
of gram-negative bacteria and this reaction is an important step in its
catalytic cycle.Stelzl  _et al._ (2020)

Figure 5: Comparison of torsion angles via the Jensen-Shannon Distance (JSD)
reveals where conformational ensembles differ and allows to visualize these
effects on a single structure. Analysis of backbone and side-chain torsions
highlights regions of nDsbD that are most affected by the oxidation of the two
cysteine thiols, Cys103 and Cys109 Stelzl  _et al._ (2020) For each residue,
the maximum JSD of the distributions of backbone torsion angles between
oxidised and reduced condition is encoded in the color of the cartoon
representation from white (0.0) via yellow and orange to red (0.3). Similarly,
the maximum JSD for the sidechain torsion angles is encoded in the color of
the side chains, where stick representation of side chains is only displayed
for residues with JSD values greater than 0.3. Structures are visualized using
PyMol.Schrödinger, LLC (2015)

Visualization of the maximum JSD per residue for backbone and sidechain
torsions (see Fig. 5) shows at one glance the opposing residues in the gap
loop to be the most affected regions. Unsurprisingly, the residues directly
involved in the reaction, Cys103 and Cys109, show the highest values. The
effects of this reaction on the surrounding residues – showing up in our
analysis as medium JSD values – cause a change in the distance of the opposing
loop and a corresponding opening of the gap. In particular, we find residues
Phe70 and Tyr71 at the neighboring loop to be strongly affected. Indeed, the
authors of the original study identified the distance between residues Phe70
and Cys109 to be the characteristic hallmark of loop opening.

Besides the influence on the gap loop, we identify a more subtly involved and
previously not discussed region at the opposite side of the protein. This
region – mainly the backbones of residues Ala92 and Gly93 – functions as a
hinge for the beta-strand region that slightly tilts when the gap loop is
pushed outward during the gap opening. This finding demonstrates how our
approach picks up differences between simulations that are otherwise easily
missed.

###  III.2 Comparing force field parameters: Interactions of Calcium with DNA

As a second use case, we show how to quantify the effects of small changes in
force field parameters on the overall conformational ensemble (Fig. 6). We
consider the binding of calcium ions (Ca2+) to DNA Cruz-León  _et al._ (2021).
Metal cations play a crucial role in stabilizing the structure of nucleic acid
systems. Their force field parameters are usually determined to reproduce bulk
properties like the solvation-free energy and thus often not directly
transferable to interactions with biomolecules.Panteva  _et al._ (2015)
Furthermore, small changes in the interactions can have significant
conformational consequences overall.Yoo and Aksimentiev (2012) Studying the
effect of such small changes on the overall conformation is an important
problem in force field optimization.

Figure 6: Comparison of torsion angles via the Kolmogorov-Smirnoff Statistic
(KSS) reveals where conformational ensembles differ and allows to visualize
these effects on a single structure. Comparison of DNA backbone torsion angles
from MD simulations to an experimental reference ensemble (see main text for
details) shows the improvement in the ensemble generated by MD simulations
when switching from standard force field parameters (left) to optimized force
field parameters (right). Structures are visualized using PyMol.Schrödinger,
LLC (2015)

Comparison of DNA backbone torsions from MD simulations using the standard
force field parameters for Ca2+-DNA interactions to an experimental reference
ensemble reveals a periodic pattern of strong deviations along the entire
double strand that are strongly reduced by using optimized force field
parameters. We generate the reference distribution from the experimental
structure by Gaussian sampling of the relevant coordinates using the
uncertainty of the X-ray structure (1.7 Å, PDB: 477D) as the width of the
distribution. We then compare the simulations of each parameter set to this
reference (Fig. 6) using the Kolmogorov-Smirnov statistic (Eq. 3) because it
provides a parameter-free measure for the deviations. The approximate
periodicity of the deviations (Fig. 6, left) shows that, using standard
parameters, the conformational ensemble as a whole deviates from the reference
distribution and has problems reproducing the overall structure of the double
strand. The authors of the original study Cruz-León  _et al._ (2021)
identified an overestimation of Ca2+-DNA interactions as the main cause of
such deviations. It allows the Ca2+ ions to bridge between the phosphate
oxygen atoms of opposite backbone strands which causes the minor groove of the
DNA strand to shrink and in turn affects the entire structure. Thus, they
rescaled the force field parameters to optimize the Ca2+-DNA interactions.
This not only improved the local accuracy but the entire conformational
ensemble, even though the overall structure was not explicitly optimized for
during the rescaling. In our PENSA-based analysis, this improvement is
immediately visible by the reduced deviations (Fig. 6, right). The overall
ensemble KSS (as in Eq. 4) is reduced from 0.50 (standard force field) to 0.34
(optimized force field). This example shows that our workflow quickly
identifies whether and where small changes in local interactions propagate to
strong deviations in the overall structure of a biomolecule.

###  III.3 Tracing information linked to a protonation state: The central
aspartic acid in the $\mu$-opioid receptor

To showcase the usefulness of State-Specific Information (SSI), we investigate
the relationship between the protonation state of a central aspartic acid and
the $\mu$-opioid receptor ($\mu$OR) ensemble. The $\mu$OR, a G protein–coupled
receptor (GPCR), is a transmembrane receptor protein that converts
extracellular stimuli into intracellular signalling cascades. The diversity in
structure and function among GPCRs underpins complex activation mechanisms
that, despite large pharmaceutical interest, remain unresolved. Katritch  _et
al._ (2013); Hauser  _et al._ (2018); Eiger  _et al._ (2022) Rotamer changes
in residue side chains give rise to larger-scale conformational changes that
enable the binding of effector proteins and trigger downstream signaling
(receptor activation).Hauser  _et al._ (2021) These rotamers can be understood
as molecular micro-switches, making them an ideal use case for our state-based
mutual information approach, SSI. Many receptors are influenced by
environmental pH changes,Ludwig  _et al._ (2003); Rowe  _et al._ (2021);
Sanderlin  _et al._ (2015) and protonation changes of an evolutionarily
conserved aspartic acid residue have been hypothesised to represent a key step
in receptor activation. Ghanouni  _et al._ (2000); Mahalingam  _et al._
(2008); Vickery  _et al._ (2018) Here, we apply SSI to side-chain and backbone
torsions to investigate the effect of protonating this conserved aspartic
acid, Asp114 ($i=$ protonated, $j=$ de-protonated) in the antagonist-bound
$\mu$OR (PDB: 4DKL, Fig. 7).

Figure 7: Analysis of backbone and side-chain torsion angles (A, cartoon and
stick), and intracavity water molecules (B, spheres), via State-Specific
Information (SSI) quantifies the information that each feature signals about
the protonation state of Asp114 in the $\mu$-opioid receptor. SSI highlights
that changes in the protonation state of Asp114 couple to rotameric state
changes on: the backbone of TM6, proximal to Asp114; the backbone of ICL2 and
ICL3; residue side chains on the P-I-F, NPxxY and DRY motifs; and water sites
distributed throughout the receptor. The SSI for backbone torsions is encoded
in the color and width of the cartoon representation from white (0.0), via
yellow (0.33) and orange (0.67), to red (1.0). Similarly, the SSI for
sidechain torsions is encoded in the color of the stick representation, and
for intracavity water sites in the color of the sphere representation. For
ease of visualization, only side chains corresponding to the twenty largest
SSI values are displayed and backbone torsions are represented solely by
cartoon. The right panel shows a close-up view of the region proximal to
Asp114, with SSI values for backbone and side chains projected onto the
receptor in the same manner as in panel A. Structures are visualized using
PyMol.Schrödinger, LLC (2015)

_Mutual information analysis of side-chain rotamers:_ The SSI values
calculated for each residue reveal those parts of the receptor that signal
information about the protonation state of Asp114 (Fig. 7). Namely, rotameric
state changes in the backbone torsions of transmembrane helix TM6, proximal to
Asp114, are coupled to the protonation state changes of Asp114. Similarly, the
backbone rotamer states of intracellular loops ICL2 and ICL3 couple to the
protonation state of Asp114. An outward swing of TM6, enabled by backbone
conformational state changes near Phe289 (P-I-F motif), is characteristic of
receptor activation,Katritch  _et al._ (2013) and specific conformations of
ICL2 and ICL3 are implicated in the binding of signal proteins.Koehl  _et al._
(2018); Mafi  _et al._ (2020) Furthermore, side chain changes are identified
on the P-I-F motif, the NPxxY motif, and the DRY motif, three receptor motifs
that are known to undergo distinct rotamer changes in the transition from
inactive to active receptor states.Katritch  _et al._ (2013) The recognition
of receptor regions where conformational changes are associated with
activation and signaling suggests that the Asp114 protonation state and GPCR
activation are intertwined. This example thus demonstrates how SSI and its
visualization help to pinpoint receptor regions where the features’ rotamer
states inform about an aspartic acids protonation state.

_Mutual information analysis of intracavity water sites:_ To further
demonstrate how SSI can be used to analyze water cavity features, we
featurized the ten most well-defined water sites (Fig. 7). The locations of
the water cavities were determined using the PENSA water featurizer as the ten
sites with the largest probability maxima in the water density grid of the
combined ensembles and labelled O1–O10 according to their ranking. The
positioning of all ten sites agrees well with water molecules resolved in
experimentsGranier  _et al._ (2012); Huang  _et al._ (2020) and predicted by
the HomolWat serverMayol  _et al._ (2020) for the inactive $\mu$OR crystal
structure (PDB: 4DKL), confirming the accuracy of PENSA’s water cavity
featurization. Three water cavities are within the vicinity of Asp114: O3, O8
and O10. Using SSI (eq. 5), we calculated that water cavities O1–O10 share
information with the Asp 114 protonation state on levels between 0.00–0.74
bits. Water cavity O8, for example, shares 0.67 bits of information in coupled
conformational state changes linked to the transition between ensemble $i$ and
$j$, i.e., unprotonated to protonated Asp114. O8 also sits beside Phe289 of
the P-I-F motif, where we identified that backbone rotamer state changes are
coupled to Asp114 protonation. Surprisingly, the more distant water site O1
shares 0.74 bits of information, 82% of which is due to an occupancy change.
Comparing the average ensemble structures reveals an increased packing between
TM5 and TM6 in the region about the water O1 cavity, with the distance between
surrounding C$\alpha$s moving over on average 1Å closer in the
Asp114-protonated ensemble, suggesting that helix movements on TM6 lead to a
collapse of the cavity. This analysis highlights a concerted behaviour of
water cavities and TM6, whereby state changes to both are indicative of the
protonation state of Asp114. It shows how the combined analysis of multiple
different features and a comprehensive visualization help to find
interrelations within a receptor and discover signaling pathways.

##  IV Discussion

With PENSA, we have implemented an open-source library that provides
systematic, easy-to-apply methods that make the otherwise often cumbersome
exploration of biomolecular systems faster, more reliable and easier to
interpret. The data readers implemented in PENSA engage with a range of
biomolecular systems via robust featurization implementations. These include
readers for interatomic distances and the characteristic torsion angles of
amino acids and nucleic acids as well as a novel approach that reliably
incorporates water and ion cavities via their occupancy and, in the case of
water, the polarization of the cavity. Combined with dimensionality reduction
tools, PENSA can handle biomolecular systems on a wide range of scales and
resolutions. PENSA includes two comparison measures with slightly different
use cases: JSD is a fast way to compare ensembles and pick up relevant
features. The sensitivity of JSD depends on — and can be adjusted via — the
spacing of the discretization. KSS is a parameter-free metric that picks up
every ensemble difference, even the noise. In most practical cases, JSD and
KSS provide similar results. In addition, PENSA includes a mutual information
measure, State-Specific information (SSI). With no prior knowledge of the
data, PENSA performs an automatic discretization of feature distributions into
conformational states, and via SSI, quantifies the information that each
features’ conformational states signal about the ensemble they are in. The
ideal use case involves features that switch between well-defined states, such
as molecular switches, however customizable state definitions allow SSI to
operate with many kinds of feature discretization. SSI can be further extended
to three or more features to quantify information flow within a system but is
currently only implemented for two. Combined with PENSA’s convenient
visualization tools, these methods allow for a detailed analysis of
biomolecular ensembles without the bias of hand-picked metrics, acting as a
solid basis for mechanistical interpretations and further, more detailed
analysis.

Our example analyses demonstrate the versatility of PENSA on three different
biomolecular systems. By investigating the effects of a small chemical
modification on loop opening in an oxidoreductase with JSD, we demonstrate the
validity of the method in confirming previously discovered results, while
additionally reporting novel, more subtle findings within the same system. We
demonstrate the applicability of PENSA in the optimization of force field
parameters via a comparison of the interactions between Calcium and DNA under
different force field parameters with KSS. Finally, we report on a
communication channel in the $\mu$-opioid receptor that transmits information
between the intracellular signalling site and the protonation state of a
distal aspartic acid, shedding further light on the signal transduction
mechanism of this mechanistically complicated system. The major limit to the
accuracy of PENSA is the quality of the input ensembles. For example,
insufficiently converged MD simulations can cause false positives when an
equally probable transition happens only in one of the conditions. Or they can
cause false negatives for overall rare events or slow processes. If in doubt,
validation by other means may be necessary (experiment, independent/longer
simulations). Although non-converged simulations can give useful hints, these
cases demand a cautious systematic analysis. Despite these caveats, PENSA has
the potential for high-throughput analysis of a large amount of simulations,
e.g., those available in GPCRmdRodríguez-Espigares  _et al._ (2020), can be
used to independently quantify the quality of force fields or generative
machine learning models, and unravel molecular mechanisms and signaling
pathways.

##  V Conclusions

In summary, we present a powerful toolkit to build workflows for the
systematic and quantitative analysis of biomolecular systems and their
conformational ensembles. PENSA provides flexible options to featurize various
biomolecular systems, metrics to compare ensembles and to detect
interrelations between different regions of a system, and methods to produce
intuitive visualizations. We demonstrate the effectiveness of these methods on
three real-world examples from molecular biology, showing how PENSA makes it
easier for researchers to analyze large amounts of complex simulation data.

##  VI Acknowledgments

We thank Lukas Stelzl for the trajectory data of the oxidoreductase example as
well as Sergio Cruz-León, Kara Grotz, and Nadine Schwierz for the trajectory
data of the DNA-Calcium example. We are grateful to Alexander Powers, Lukas
Stelzl, Nicole Ong, Eleanore Ocana, Emma Andrick, Callum Ives, and Bu Tran for
beta-testing as well as Maria Karelina, Marc Dämgen, Patricia Suriana, Sergio
Cruz-León, Michael Ward, and Ramon Guixà-González for helpful discussions.
M.V. was supported by the EMBO long-term fellowship ALTF 235-2019. N.J.T. was
supported by a BBSRC EASTBIO PhD studentship. This work was supported by
National Institutes of Health grant R01GM127359 (R.O.D.) An award of computer
time was provided by the INCITE program. This research used resources of the
Oak Ridge Leadership Computing Facility, a DOE Office of Science User Facility
supported under contract DE-AC05-00OR22725. Additional computing for this
project was performed on the Sherlock cluster and the University of Dundee SLS
HPC cluster. We thank Stanford University, the Stanford Research Computing
Facility, and the University of Dundee for providing computational resources
and support that contributed to these research results.

## References

  * Hollingsworth and Dror (2018) S. A. Hollingsworth and R. O. Dror, Neuron 99, 1129 (2018).
  * Knoverek  _et al._ (2018) C. R. Knoverek, G. K. Amarasinghe, and G. R. Bowman, Trends in Biochemical Sciences 44, 351 (2018).
  * Suomivuori  _et al._ (2020) C.-M. Suomivuori, N. R. Latorraca, L. M. Wingler, S. Eismann, M. C. King, A. L. W. Kleinhenz, M. A. Skiba, D. P. Staus, A. C. Kruse, R. J. Lefkowitz, and R. O. Dror, Science 367, 881 (2020).
  * Provasi  _et al._ (2011) D. Provasi, M. C. Artacho, A. Negri, J. C. Mobarec, and M. Filizola, PLoS Computational Biology 7, 1 (2011).
  * McCorvy  _et al._ (2018) J. D. McCorvy, K. V. Butler, B. Kelly, K. Rechsteiner, J. Karpiak, R. M. Betz, B. L. Kormos, B. K. Shoichet, R. O. Dror, J. Jin, and B. L. Roth, Nature Chemical Biology 14, 126 (2018).
  * Cordero-Morales  _et al._ (2007) J. F. Cordero-Morales, V. Jogini, A. Lewis, V. Vásquez, D. M. Cortes, B. Roux, and E. Perozo, Nature Structural and Molecular Biology 14, 1062 (2007).
  * Liu  _et al._ (2015) Y. Liu, M. Ke, and H. Gong, Biophysical Journal 109, 542 (2015).
  * Thomson  _et al._ (2021) N. J. Thomson, O. N. Vickery, C. M. Ives, and U. Zachariae, bioRxiv (2021), 10.1101/2020.08.28.271510.
  * Humphrey  _et al._ (1996) W. Humphrey, A. Dalke, and K. Schulten, Journal of Molecular Graphics 14, 33 (1996).
  * DeLano (2002) W. L. DeLano, CCP4 Newsletter on protein crystallography 40, 82 (2002).
  * Dror  _et al._ (2011) R. O. Dror, D. H. Arlow, P. Maragakis, T. J. Mildorf, A. C. Pan, H. Xu, D. W. Borhani, and D. E. Shaw, Proceedings of the National Academy of Sciences of the United States of America 108, 18684 (2011).
  * Zhou  _et al._ (2019) Q. Zhou, D.-H. Yang, M. Wu, Y. Guo, W. Guo, L. Zhong, X. Cai, A. Dai, W. Jang, E. I. Shakhnovich, Z.-J. Liu, R. C. Stevens, N. A. Lambert, M. M. Babu, M.-W. Wang, and S. Zhao, eLife 8, e50279 (2019).
  * Dror  _et al._ (2013) R. O. Dror, H. F. Green, C. Valant, D. W. Borhani, J. R. Valcourt, A. C. Pan, D. H. Arlow, M. Canals, J. R. Lane, R. Rahmani, J. B. Baell, P. M. Sexton, A. Christopoulos, and D. E. Shaw, Nature 503, 295 (2013).
  * Bowman  _et al._ (2015) G. R. Bowman, E. R. Bolin, K. M. Hart, B. C. Maguire, and S. Marqusee, Proceedings of the National Academy of Sciences of the United States of America 112, 2734 (2015).
  * Zivanovic  _et al._ (2020) S. Zivanovic, G. Bayarri, F. Colizzi, D. Moreno, J. L. Gelpí, R. Soliva, A. Hospital, and M. Orozco, Journal of Chemical Theory and Computation 16, 6586 (2020).
  * Brüschweiler (2003) R. Brüschweiler, Proteins: Structure, Function and Genetics 50, 26 (2003).
  * Lindorff-Larsen and Ferkinghoff-Borg (2009) K. Lindorff-Larsen and J. Ferkinghoff-Borg, PLOS ONE 4, 1 (2009).
  * Tiberti  _et al._ (2015) M. Tiberti, E. Papaleo, T. Bengtsen, W. Boomsma, and K. Lindorff-Larsen, PLoS Computational Biology 11, 1 (2015).
  * Noé  _et al._ (2019) F. Noé, S. Olsson, J. Köhler, and H. Wu, Science 365, eaaw1147 (2019).
  * Husic and Pande (2018) B. E. Husic and V. S. Pande, Journal of the American Chemical Society 140, 2386 (2018).
  * Nüske  _et al._ (2017) F. Nüske, H. Wu, J. H. Prinz, C. Wehmeyer, C. Clementi, and F. Noé, Journal of Chemical Physics 146, 094104 (2017).
  * Fraccalvieri  _et al._ (2011) D. Fraccalvieri, A. Pandini, F. Stella, and L. Bonati, BMC Bioinformatics 12, 158 (2011).
  * Ward  _et al._ (2021) M. D. Ward, M. I. Zimmerman, A. Meller, M. Chung, S. J. Swamidass, and G. R. Bowman, Nature Communications 12, 1 (2021).
  * Malik  _et al._ (2021) M. Malik, M. D. Ward, Y. Fang, J. R. Porter, M. I. Zimmerman, T. Koelblen, M. Roh, A. I. Frolova, T. P. Burris, G. R. Bowman, P. I. Imoukhuede, and S. K. England, ACS Pharmacology and Translational Science 4, 1543 (2021).
  * Vögele  _et al._ (2021) M. Vögele, N. Thomson, S. Truong, and J. McAvity, “PENSA,”  (2021), doi: 10.5281/ZENODO.4362136.
  * Schrödinger, LLC (2015) Schrödinger, LLC, “The PyMOL molecular graphics system, version 1.8,” (2015).
  * Kunzmann and Hamacher (2018) P. Kunzmann and K. Hamacher, BMC Bioinformatics 19, 346 (2018).
  * Stelzl  _et al._ (2020) L. S. Stelzl, D. A. Mavridou, E. Saridakis, D. Gonzalez, A. J. Baldwin, S. J. Ferguson, M. S. Sansom, and C. Redfield, eLife 9, 1 (2020).
  * Cruz-León  _et al._ (2021) S. Cruz-León, K. K. Grotz, and N. Schwierz, Journal of Chemical Physics 154, 171102 (2021).
  * Michaud-Agrawal  _et al._ (2011) N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein, Journal of Computational Chemistry 32, 2319 (2011).
  * Gowers  _et al._ (2016) R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, J. Domański, D. L. Dotson, S. Buchoux, I. M. Kenney, and O. Beckstein, in  _Proceedings of the 15th Python in Science Conference_ (2016) pp. 98–105.
  * Scherer  _et al._ (2015) M. K. Scherer, B. Trendelkamp-Schroer, F. Paul, G. Pérez-Hernández, M. Hoffmann, N. Plattner, C. Wehmeyer, J. H. Prinz, and F. Noé, Journal of Chemical Theory and Computation 11, 5525 (2015).
  * Keating  _et al._ (2011) K. S. Keating, E. L. Humphris, and A. M. Pyle, Quarterly Reviews of Biophysics 44, 433 (2011).
  * Venkatakrishnan  _et al._ (2019) A. J. Venkatakrishnan, A. K. Ma, R. Fonseca, N. R. Latorraca, B. Kelly, R. M. Betz, C. Asawa, B. K. Kobilka, and R. O. Dror, Proceedings of the National Academy of Sciences of the United States of America 116, 3288 (2019).
  * Yuan  _et al._ (2014) S. Yuan, S. Filipek, K. Palczewski, and H. Vogel, Nature Communications 5, 4733 (2014).
  * Pardo  _et al._ (2007) L. Pardo, X. Deupi, N. Dölker, M. L. López-Rodríguez, and M. Campillo, ChemBioChem 8, 19 (2007).
  * Levy and Onuchic (2004) Y. Levy and J. N. Onuchic, Proceedings of the National Academy of Sciences of the United States of America 101, 3325 (2004).
  * Kajander  _et al._ (2000) T. Kajander, P. C. Kahn, S. H. Passila, D. C. Cohen, L. Lehtiö, W. Adolfsen, J. Warwicker, U. Schell, and A. Goldman, Structure 8, 1203 (2000).
  * Zarzycka  _et al._ (2019) B. Zarzycka, S. A. Zaidi, B. L. Roth, and V. Katritch, Pharmacological Reviews 71, 571 (2019).
  * Andreini  _et al._ (2008) C. Andreini, I. Bertini, G. Cavallaro, G. L. Holliday, and J. M. Thornton, Journal of Biological Inorganic Chemistry 13, 1205 (2008).
  * Ives  _et al._ (2022) C. M. Ives, N. J. Thomson, and U. Zachariae, bioRxiv (2022), 10.1101/2022.04.01.486690.
  * Dunbrack (2002) R. L. Dunbrack, Current Opinion in Structural Biology 12, 431 (2002).
  * Dunbrack and Karplus (1993) R. L. Dunbrack and M. Karplus, Journal of Molecular Biology 230, 543 (1993).
  * Scouras and Daggett (2011) A. D. Scouras and V. Daggett, Protein Science 20, 341 (2011).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * F.R.S. (1901) K. P. F.R.S., The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2, 559 (1901).
  * Schultze and Grubmüller (2021) S. Schultze and H. Grubmüller, Journal of Chemical Theory and Computation 17, 5766 (2021).
  * Macqueen (1967) J. Macqueen, in  _In 5-th Berkeley Symposium on Mathematical Statistics and Probability_ (1967) pp. 281–297.
  * Prinz  _et al._ (2011) J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D. Chodera, C. Schütte, and F. Noé, The Journal of Chemical Physics 134, 174105 (2011).
  * Shannon (1948) C. E. Shannon, Bell System Technical Journal 27, 379 (1948).
  * Bell (2003) A. J. Bell, 4th International Symposium on Independent Component Analysis and Blind Source Separation (2003).
  * Hunter (2007) J. D. Hunter, Computing in Science and Engineering 9, 90 (2007).
  * Panteva  _et al._ (2015) M. T. Panteva, G. M. Giambaşu, and D. M. York, Journal of Physical Chemistry B 119, 15460 (2015).
  * Yoo and Aksimentiev (2012) J. Yoo and A. Aksimentiev, Journal of Physical Chemistry Letters 3, 45 (2012).
  * Katritch  _et al._ (2013) V. Katritch, V. Cherezov, and R. C. Stevens, Annual Review of Pharmacology and Toxicology 53, 531 (2013).
  * Hauser  _et al._ (2018) A. S. Hauser, S. Chavali, I. Masuho, L. J. Jahn, K. A. Martemyanov, D. E. Gloriam, and M. M. Babu, Cell 172, 41 (2018).
  * Eiger  _et al._ (2022) D. S. Eiger, U. Pham, J. Gardner, C. Hicks, and S. Rajagopal, American Journal of Physiology-Cell Physiology 322, C887 (2022).
  * Hauser  _et al._ (2021) A. S. Hauser, A. J. Kooistra, C. Munk, F. M. Heydenreich, D. B. Veprintsev, M. Bouvier, M. M. Babu, and D. E. Gloriam, Nature Structural and Molecular Biology 28, 879 (2021).
  * Ludwig  _et al._ (2003) M. G. Ludwig, M. Vanek, D. Guerini, J. A. Gasser, C. E. Jones, U. Junker, H. Hofstetter, R. M. Wolf, and K. Seuwen, Nature 425, 93 (2003).
  * Rowe  _et al._ (2021) J. B. Rowe, N. J. Kapolka, G. J. Taghon, W. M. Morgan, and D. G. Isom, Journal of Biological Chemistry 296, 100167 (2021).
  * Sanderlin  _et al._ (2015) E. J. Sanderlin, C. R. Justus, E. A. Krewson, and L. V. Yang, Cell Health and Cytoskeleton 7, 99 (2015).
  * Ghanouni  _et al._ (2000) P. Ghanouni, H. Schambye, R. Seifert, T. W. Lee, S. G. Rasmussen, U. Gether, and B. K. Kobilka, Journal of Biological Chemistry 275, 3121 (2000).
  * Mahalingam  _et al._ (2008) M. Mahalingam, K. Martínez-Mayorga, M. F. Brown, and R. Vogel, Proceedings of the National Academy of Sciences of the United States of America 105, 17795 (2008).
  * Vickery  _et al._ (2018) O. N. Vickery, C. A. Carvalheda, S. A. Zaidi, A. V. Pisliakov, V. Katritch, and U. Zachariae, Structure 26, 171 (2018).
  * Koehl  _et al._ (2018) A. Koehl, H. Hu, S. Maeda, Y. Zhang, Q. Qu, J. M. Paggi, N. R. Latorraca, D. Hilger, R. Dawson, H. Matile, G. F. X. Schertler, S. Granier, W. I. Weis, R. O. Dror, A. Manglik, G. Skiniotis, B. K. Kobilka, C. J. Draper-Joyce, M. Khoshouei, D. M. Thal, Y.-L. Liang, A. T. N. Nguyen, S. G. B. Furness, H. Venugopal, J.-A. Baltos, J. M. Plitzko, R. Danev, and W. Baumeister, Nature , 547–552 (2018).
  * Mafi  _et al._ (2020) A. Mafi, S. K. Kim, and W. A. Goddard, Proceedings of the National Academy of Sciences of the United States of America 117, 16346 (2020).
  * Granier  _et al._ (2012) S. Granier, A. Manglik, A. C. Kruse, T. S. Kobilka, F. S. Thian, W. I. Weis, B. K. Kobilka, J. M. Mathiesen, R. K. Sunahara, and L. Pardo, Nature 485, 321–326 (2012).
  * Huang  _et al._ (2020) W. Huang, A. Manglik, A. J. Venkatakrishnan, T. Laeremans, E. N. Feinberg, A. L. Sanborn, H. E. Kato, K. E. Livingston, T. S. Thorsen, R. C. Kling, S. Granier, P. Gmeiner, S. M. Husbands, J. R. Traynor, W. I. Weis, J. Steyaert, R. O. Dror, and B. K. Kobilka,  _Nature_ , Nature 584, E16 (2020).
  * Mayol  _et al._ (2020) E. Mayol, A. Garcia-Recio, J. K. Tiemann, P. W. Hildebrand, R. Guixa-Gonzalez, M. Olivella, and A. Cordomi, Nucleic Acids Research 48, W54 (2020).
  * Rodríguez-Espigares  _et al._ (2020) I. Rodríguez-Espigares, M. Torrens-Fontanals, J. K. Tiemann, D. Aranda-García, J. M. Ramírez-Anguita, T. M. Stepniewski, N. Worp, A. Varela-Rial, A. Morales-Pastor, B. Medel-Lacruz, G. Pándy-Szekeres, E. Mayol, T. Giorgino, J. Carlsson, X. Deupi, S. Filipek, M. Filizola, J. C. Gómez-Tamayo, A. Gonzalez, H. Gutiérrez-de Terán, M. Jiménez-Rosés, W. Jespers, J. Kapla, G. Khelashvili, P. Kolb, D. Latek, M. Marti-Solano, P. Matricon, M. T. Matsoukas, P. Miszta, M. Olivella, L. Perez-Benito, D. Provasi, S. Ríos, I. R. Torrecillas, J. Sallander, A. Sztyler, S. Vasile, H. Weinstein, U. Zachariae, P. W. Hildebrand, G. De Fabritiis, F. Sanz, D. E. Gloriam, A. Cordomi, R. Guixà-González, and J. Selent, Nature Methods 17, 777 (2020).


# An Empirical Study of Protocols in Smart Contracts

Timothy Mou<EMAIL_ADDRESS>Swarthmore CollegeSwarthmorePennsylvaniaUSA
,  Michael Coblenz<EMAIL_ADDRESS>University of MarylandCollege
ParkMarylandUSA and  Jonathan Aldrich<EMAIL_ADDRESS>Carnegie
Mellon UniversityPittsburghPennsylvaniaUSA

(2021)

###### Abstract.

Smart contracts are programs that are executed on a blockhain. They have been
used for applications in voting, decentralized finance, and supply chain
management. However, vulnerabilities in smart contracts have been abused by
hackers, leading to financial losses. Understanding state machine protocols in
smart contracts has been identified as important to catching common bugs,
improving documentation, and optimizing smart contracts. We analyze Solidity
smart contracts deployed on the Ethereum blockchain and study the prevalence
of protocols and protocol-based bugs, as well as opportunities for gas
optimizations.

smart contracts, Ethereum, Solidity

††copyright: acmcopyright††journalyear: 2021††doi:
10.1145/1122445.1122456††booktitle: HATRA ’21: Human Aspects of Types and
Reasoning Assistants, October 19, 2021, Chicago, IL††ccs: Software and its
engineering Distributed programming languages††ccs: Software and its
engineering State systems

##  1\. Introduction

Smart contracts are general-purpose programs that are run on a blockchain and
maintain state stored on a blockchain’s ledger. They have been used in a
variety of applications, including voting, decentralized finance, and supply
chain management (IBM, 2019; Elsden et al., 2018). However, vulnerabilities in
smart contracts can be exploited and have serious consequences, such as the
well-known DAO hack in which over $50 million worth of cryptocurrency was
stolen (Popper, 2016). Smart contracts on Ethereum cannot be modified once
they are deployed, making bugs difficult to fix retroactively. In addition, on
Ethereum, the most common platform for smart contracts, users must pay gas in
order to run transactions, which is used to reward miners and prevent non-
terminating transactions. High gas prices can make interacting with smart
contracts prohibitively expensive, limiting their usefulness (Rozen, 2021).

Smart contracts have been observed to often follow protocols that can be
represented by state machines (Foundation, 2021). Protocols specify that a
contract’s functions must be called in a particular order. A contract may have
multiple abstract states that it transitions between, where each abstract
state determines which functions can be called in that particular state. For
instance, consider a wallet contract whose public functions are shown in
Figure 1. Its purpose is to transfer its token balance to a beneficiary after
a fixed amount of time. The wallet begins in an initial state, and to initiate
the transfer, the owner of the wallet must call lock() to lock the wallet for
a set period of time. Only after lock() has been called and the time has
elapsed can this state can release() be called to transfer the wallet’s
balance to the beneficiary.

⬇

1function lock() public onlyOwner returns (bool) {

2 require(!isLocked);

3 require(tokenBalance() > 0);

4 start_time = now;

5 end_time = start_time.add(fifty_two_weeks);

6 isLocked = true;

7}

8

9function lockOver() constant public returns (bool) {

10 uint256 current_time = now;

11 return current_time > end_time;

12}

13

14function release() onlyOwner public {

15 require(isLocked);

16 require(!isReleased);

17 require(lockOver());

18 uint256 token_amount = tokenBalance();

19 token_reward.transfer(beneficiary, token_amount);

20 emit TokenReleased(beneficiary, token_amount);

21 isReleased = true;

22}

Figure 1. The main functions in a wallet contract. The onlyOwner modifier on
the lock and release functions allows only the owner of the wallet contract to
call them.

The concept of typestate (Strom and Yemini, 1986) allows one to specify the
set of valid operations which can be called on a contract, which may change as
the contract transitions to different states. New programming languages for
smart contracts, including Obsidian (Coblenz et al., 2020) and Flint (Schrans
et al., 2018), have proposed that typestate can be used to express protocols
more explicitly and statically check adherence to protocols. One benefit of
such systems is that they allow developers to specify functions and fields
that are only in scope in certain states. Failure to change or reset fields
between state transitions can lead to bugs that are difficult to track down.
In addition, accessing storage on the blockchain is costly, and using
typestate may allow contracts to lower their gas usage by clearing fields or
avoiding unnecessary writes to storage.

We conducted an empirical study of Solidity smart contracts to demonstrate
that protocols are commonly used in smart contracts and that typestate can
help contracts avoid bugs and use less gas. We collected a sample of 100 smart
contracts and determined that 37 of them defined protocols. We wrote a static
analysis tool, StateOptimizationDetector, to find opportunities to lower the
cost of transactions by taking advantage of Ethereum’s gas refund policy. We
ran StateOptimizationDetector on a sample of 1,000 contracts and found that 80
(0.8%) of them could use this optimization, with a false positive rate of
28.6%.

We investigated the prevalence of one access control bug, which allows
previous owners of a contract to unexpectedly reclaim ownership, to show that
the incorrect implementation of protocols can create vulnerabilities. We wrote
a static analysis tool, AccessLockDetector, to identify contracts which could
have bugs in the implementation of their access control mechanisms. We ran
AccessLockDetector on a sample of 4,500 contracts and identified 36 (0.8%) of
them via manual examination as vulnerable to this exploit, with a false
positive rate of 34.5%. We observed a high level of similarity among the
identified vulnerable contracts, suggesting that many of them originated from
a common source. We hypothesize that language features that allow protocols to
be expressed more directly can help contract writers avoid these bugs or catch
them before deployment.

##  2\. Methodology

To inform our understanding of how smart contracts define and use protocols,
we first tried to identify contracts with protocols by looking for patterns in
contracts’ source code. We built a detector using the Slither static analysis
tool (Feist et al., 2019), named ProtocolDetector, that identifies contracts
with protocols by searching for reverting execution paths in functions that
resulted from reading the contract’s fields. The motivation behind this is
that when a function in a protocol-defining contract is called, the contract
will often first check that it is in an appropriate abstract state, which is
typically represented by one or more fields. Reverting paths that result from
reading fields often indicate that the protocol has been violated. Similar
patterns have been used to search for evidence of object protocols in Java
(Beckman et al., 2011). For example, in Figure 1, the function release will
revert if the values of the fields isLocked and isReleased are not true and
false, respectively. This indicates that these fields represent the contract’s
abstract state, and that this contract defines a protocol.

We incorporated several heuristics specific to Solidity to help the detector
find common patterns indicative of protocols. For example, Solidity allows the
definition of modifiers, which can be attached to function declarations to
insert code that is run before or after the main body of the function. When
checking a function for state patterns, the detector will also check the body
of any of its modifiers, since they are commonly used to ensure that the
contract is in a valid state before executing the function.

###  2.1. Finding protocol-based optimizations

We then identified a method to lower the gas cost of some transactions using
Ethereum’s refund mechanism. When a storage value is reset (set to zero from
non-zero), a refund of (up to) 15,000 gas is given back to the sender of the
transaction (Wood et al., 2014). Contracts with certain types of protocols may
be able to use the refund to lower the cost of some of their transactions.
When a contract transitions to a new state, some of its fields may no longer
be used for the remaining lifetime of the contract. These fields could
therefore be reset during this state transition to take advantage of the
refund, lowering the overall gas cost. This optimization encourages the
clearing of unused fields, reducing the amount of unneeded state that is kept
on the Ethereum blockchain (Moosavi and Clark, 2021). See Figure 2 for an
example case of this optimization. Here, the function addToX is “deactivated”
by calling stop(). Since the field amountToAdd is only used in the function
addToX, after stop() is called, amountToAdd can safely be reset without
affecting any observable behavior of the contract. It should be noted that a
recent proposal to Ethereum that is included in the London hard fork, EIP-3529
(Buterin and Swende, 2021), reduces the amount of gas that is refunded by
clearing storage, making this optimization no longer viable after the proposal
is implemented. However, this optimization is still applicable to earlier
versions of Ethereum.

⬇

1contract Deactivation {

2 bool public stopped;

3 uint public x;

4 uint private amountToAdd;

5

6 function stop() public {

7 require(!stopped);

8 stopped = true;

9 amountToAdd = 0; //This field can safely be reset here.

10 }

11

12 function addToX() public {

13 require(!stopped);

14 amountToAdd++;

15 x += amountToAdd;

16 }

17}

Figure 2. An example of a protocol-based optimization.

To investigate how useful this optimization was, we identified two categories
of protocols (described in Section 2) that indicated that this optimization
was applicable. Then we built a detector, named StateOptimizationDetector, to
find contracts with these protocols. This detector infers the contract’s
transition graph between different abstract states by identifying the set of
states a contract can be in before and directly after calling a function. The
detector uses the state transition graph to identify fields that will never be
used after a state transition, meaning that it will never be read from or
written to in any function that can be called in the current state or any
possible future state.

We manually examined each candidate contract reported by the
StateOptimizationDetector. During this process, we checked that each candidate
contract had one of the two protocol patterns, and that each identified field
could safely be reset without affecting the contract’s expected behavior. One
issue we ran into was that some fields are declared public, meaning that a
public getter function for the field is automatically generated and added to
the contract’s interface. Thus, resetting this field could affect the
contract’s clients that access this getter function after the state
transition, even if it does not affect the behavior of the contract itself.
However, the fact that this field is no longer used may also indicate that
this field is conceptually no longer in scope during this state. Whether to
count such a field depends on the specifics of the contract, and when manually
examining the results, we only counted the fields that we deemed safe to
reset. For example, in some contracts, we noted that fields would not be used
within the contract, but were simply used to maintain some status that would
be useful to other contracts. An example is shown in Figure 3. Such contracts
were discarded during the manual review process.

⬇

1contract MintableToken is StandardToken, Ownable {

2 uint public totalSupply;

3 bool public mintingFinished = false;

4

5 modifier canMint() {

6 require(!mintingFinished);

7 _;

8 }

9

10 function mint(address _to, uint256 _amount) onlyOwner canMint returns
(bool) {

11 totalSupply = totalSupply.add(_amount);

12 balances[_to] = balances[_to].add(_amount);

13 return true;

14 }

15

16 function finishMinting() onlyOwner returns (bool) {

17 mintingFinished = true;

18 return true;

19 }

20}

Figure 3. A token contract with a field totalSupply that maintains the total
number of tokens. After finishMinting() is called, the value of totalSupply
can no longer change, but it should not be reset, since the total number of
tokens should still be available even after no more tokens can be minted.

###  2.2. Finding bugs in access control contracts

⬇

1contract Ownable is Context {

2 address private _owner;

3 address private _previousOwner;

4 uint256 private _lockTime;

5

6 ...

7

8 /**

9 * @dev Throws if called by any account other than the owner.

10 */

11 modifier onlyOwner() {

12 require(_owner == _msgSender(), "Ownable: caller is not the owner");

13 _;

14 }

15

16 /**

17 * @dev Leaves the contract without owner. It will not be possible to call

18 * ‘onlyOwner‘ functions anymore. Can only be called by the current owner.

19 *

20 * NOTE: Renouncing ownership will leave the contract without an owner,

21 * thereby removing any functionality that is only available to the owner.

22 */

23 function renounceOwnership() public virtual onlyOwner {

24 emit OwnershipTransferred(_owner, address(0));

25 _owner = address(0);

26 }

27

28 /**

29 * @dev Transfers ownership of the contract to a new account (‘newOwner‘).

30 * Can only be called by the current owner.

31 */

32 function transferOwnership(address newOwner) public virtual onlyOwner {

33 require(newOwner != address(0), "Ownable: new owner is the zero address");

34 emit OwnershipTransferred(_owner, newOwner);

35 _owner = newOwner;

36 }

37

38 function geUnlockTime() public view returns (uint256) {

39 return _lockTime;

40 }

41

42 //Locks the contract for owner for the amount of time provided

43 function lock(uint256 time) public virtual onlyOwner {

44 _previousOwner = _owner;

45 _owner = address(0);

46 _lockTime = now + time;

47 emit OwnershipTransferred(_owner, address(0));

48 }

49

50 //Unlocks the contract for owner when _lockTime is exceeds

51 function unlock() public virtual {

52 require(_previousOwner == msg.sender, "You don’t have permission to
unlock");

53 require(now > _lockTime , "Contract is locked until 7 days");

54 emit OwnershipTransferred(_owner, _previousOwner);

55 _owner = _previousOwner;

56 }

57}

Figure 4. A Ownable contract with a temporary lock mechanism (lines 42-56).
Note that the unlock function allows ownership to be reclaimed after it has
been renounced or transferred to a new owner.

Access control mechanisms are used in smart contracts to restrict access to
certain administrative functions. An example of an access control contract is
the Ownable contract shown in Fig. 4, which defines an onlyOwner modifier that
can be attached to functions, which allows only the owner of a contract to
call them. Contracts can inherit from the Ownable contract and then use the
onlyOwner modifier to restrict access to certain functions.

Access control contracts also commonly have the ability to transfer ownership
to another address and renounce ownership (preventing anyone from calling
onlyOwner functions). Some access control contracts also have a temporary lock
feature, which allows the owner to temporarily relinquish ownership,
preventing themselves and anyone else from calling onlyOwner functions for a
set period of time. During this time, users of a contract can be assured that
onlyOwner functions cannot be called. A contract with both a temporary lock
and the ability to renounce ownership can be modeled with a protocol with
three abstract states, Locked, Unlocked, and Renounced.

However, we observed that many implementations of this temporary lock feature
violate this protocol and create a vulnerability in the expected behavior of
renouncing and transferring ownership. As shown in Table 1, after the function
unlock is called once, it can continue to be called even after the owner has
renounced or transferred ownership, since the value of _previousOwner is never
reset. This interaction may affect users of a contract who make decisions
based on the assumption that onlyOwner functions cannot be called after a
contract is renounced. For instance, a malicious owner of a token contract
could renounce ownership to convince prosepective buyers that the price of the
token is fixed, then use this exploit to reclaim ownership and raise the price
of the token.

Table 1. Steps to exploit a vulnerability in the Ownable contract in Fig. 4, where the owner of a contract can take back ownership after they have already renounced ownership. | _owner | _previousOwner | Contract state
---|---|---|---
Contract is initialized by an address addr. _owner is set to addr. | addr | 0x0 | Unlocked
Owner calls lock(0). | 0x0 | addr | Locked
Owner calls unlock(). | addr | addr | Unlocked
Owner calls renounceOwnership(). | 0x0 | addr | Renounced
Owner calls unlock(). | addr | addr | Unlocked

To investigate how commonly this bug occurs, we built a detector,
AccessLockDetector, which is designed to identify access control contracts
which define temporary lock mechanisms. To find access control contracts, it
looks for contracts that have an “owner” address field, which represents the
only address that can call restricted functions, and define an onlyOwner
modifier, which modifies functions so that they can only be called by the
contract’s owner. Then, it looks for a “locking” function in which the owner
address is first assigned to another field and then reset. This pattern
indicates that owner’s access has been disabled, and that the owner address
has been saved to another field so access can be re-enabled later. We observed
that contracts that implemented this temporary lock mechanism frequently
contain a bug that breaks the expected behavior of renouncing and transferring
ownership.

##  3\. Results

###  3.1. Prevalence of protocol-defining contracts

We ran the ProtocolDetector on a random sample of 100 contracts from the
SmartBugs Wild dataset (Durieux et al., 2020), a collection of over 47,398
unique Solidity contracts. ProtocolDetector identified 71 protocols in 37
contracts. We manually reviewed each reported protocol and categorized them
based on each protocol’s intended use and method of transitioning between
states. We determined that 2 (2.8%) of the protocols found were false
positives. However, in both cases, the contract defined other protocols, so in
total, 37 of the contracts defined protocols. The categorization of protocols,
which is based of a categorization of Java protocols in (Beckman et al.,
2011), is shown in Table 2.

Table 2. Manual categorization of 69 protocols identified with the ProtocolDetector from a random sample of 100 contracts. Category | Count |  Description
---|---|---
Toggle | 21 |  The owner of a contract can enable/disable functionality at will.
Activation | 12 |  Some functions become eventually enabled.
Redundant Operation | 12 |  Prevents a function from being called twice.
Deactivation | 10 |  Some functions become eventually disabled.
Threshold | 6 |  A function is enabled once a counter reaches a certain value.
Interval | 3 |  A function is enabled during a specified interval of time.
Lifecycle | 3 |  The contract transitions through various states modeling the lifetime of the contract.
Other | 2 |

###  3.2. Protocol-based optimizations

From the categories of protocols that we had observed, two of them indicated
contracts which could possibly be optimized by resetting fields. Deactivation
protocols have two states that are represented by a boolean field and have a
set of functions that are initially enabled. Once the boolean value changes,
the functions are permanently disabled. Lifecycle protocols typically have a
set of states through which the contract monotonically transitions through,
eventually ending in some terminal state. Contracts which have either
Deactivation or Lifecycle protocols permanently lose functionality as the
contract transitions, which may lead to fields which are no longer used and
can be safely reset.

We ran the StateOptimizationDetector on a random sample of 1,000 contracts
from the SmartBugs Wild dataset. The StateOptimizationDetector reported 112 of
these as candidates for optimization, 80 of which we determined were true
positives.

Table 3. Results of running the StateOptimizationDetector on a random sample of 1,000 contracts. Category | Candidates | Optimizable
---|---|---
Deactivation | 96 | 73
Lifecycle | 16 | 7
Total | 112 (11.2%) | 80 (8.0%)

###  3.3. Protocol-based access control bugs

We ran the AccessLockDetector on a random sample of 4,500 contracts from the
Smart Contract Sanctuary (Ortner and Eskandari, 2021), a collection of over
99,582 smart contracts deployed on the main Ethereum blockchain, of which the
AccessLockDetector flagged 55 of these contracts. We used the Smart Contract
Sanctuary dataset because it is updated daily, and we observed that the access
control bug we were looking for only appeared in contracts that were deployed
within the past year. We chose a sample of 4,500 contracts so that all of the
contracts flagged by the detector could be reviewed by hand. We then manually
examined each candidate contract, checking if it was an access control
contract that defined a temporary lock as well as functions for renouncing
ownership. Then we checked if its implementation allowed the unlock function
to be called when the contract was already unlocked, creating the
vulnerability where ownership could be reclaimed after renouncing ownership.
After examining the positive results, we determined that 36 of them were
access control contracts that defined a temporary lock feature and had this
vulnerability.

All of the 36 true positives followed the same structure as in Figure 4 in
defining the lock and unlock functions. 34 of these contracts contained a
misleading error message in the unlock function (line 53 of Figure 4), and 33
of them contained the same typo in a function name (line 38 of Figure 4). In
addition, all of the contracts passed an incorrect argument to the
OwnershipTransferred event in the lock function (line 47 of Figure 4; the
first argument should be _previousOwner, since _owner has already been reset).
Although events are not used by other smart contracts, this could provide
incorrect information to applications outside the blockchain that are
subscribed to this event.

##  4\. Related Work

Previous studies have employed static analysis and symbolic execution tools to
analyze bugs and security vulnerabilities in smart contracts. We built our
detectors using Slither (Feist et al., 2019), an extensible static analysis
tool which analyses Solidity source code and converts it into a convenient
intermediate representation. Oyente (Luu et al., 2016) is a symbolic execution
tool that finds common bugs involving issues with transaction-ordering
dependence and reentrancy attacks. Smartcheck (Tikhomirov et al., 2018) runs a
syntactical analysis on Solidity source code to detect a number of code issues
and poor practices.

A high level of code cloning and duplication among smart contracts has been
shown through various studies and smart contract datasets. In the SmartBugs
Wild dataset (Durieux et al., 2020), only 47,398 out of 972,855 (4.8%)
verified contracts (those with source code available) were unique. In a
dataset of 22,493 smart contracts used for (Kalra et al., 2018), only 1,524
(6.8%) of contracts were unique. A study on code cloning in smart contracts
found that only 20.8% of verified contracts were not a clone of any other
verified contract (Kondo et al., 2020).

Gasper (Chen et al., 2017) is a tool that locates gas-costly patterns by
analyzing contracts’ bytecodes. It reduces the cost of transactions by
eliminating dead code and rewriting expensive loop operations. We are not
aware of any tools which optimize gas usage by taking advantage of Ethereum’s
gas refund policy.

##  5\. Discussion and Future Work

We have shown that protocols appear relatively frequently in Solidity smart
contracts, with 37% of contracts defining protocols. For comparison, it is
estimated that 7.2% of types in Java programs define protocols (Beckman et
al., 2011). Therefore, we believe the implementation of protocols is an
important consideration for the design of languages for smart contracts.
Moreover, a number of categories of protocols allow for gas optimizations,
where fields can safely be reset to reduce the cost of transactions.

We have demonstrated that approximately 0.8% of deployed Ethereum contracts
contain a bug in their access control mechanisms that stems from a failure to
accurately model a contract’s protocol. We believe that introducing language
features such as typestate could help developers avoid writing these bugs like
this in the first place. Analyzing the Ownable contract from Figure 4, we
notice that the vulnerability arises from the fact that the unlock function
can be called even when the contract is not locked because the value of
_previousOwner is never reset. In a language such as Obsidian, one can
describe a contract as having multiple abstract states, and functions and
fields can be specified as being in scope only in certain states. This allows
the expected behavior of the Ownable to be expressed more easily: the Ownable
contract has three states, Locked, Unlocked, and Renounced, and the unlock
function is only in scope in the Locked state. In addition, one can avoid
having to define a _previousOwner field altogether simply by specifying the
onlyOwner modifier to only be in scope in the Unlocked state, preventing
onlyOwner functions from being called if the contract is in the Locked or
Renounced states.

We found that many of the vulnerable contracts flagged by the
AccessLockDetector had nearly identical code in their access control functions
and contained the same mistakes in function names and error messages. This
supports previous results about the high level of code cloning and duplication
between smart contracts. We hypothesize that these mistakes indicate that
these contracts were copied from a single source, rather than created
independently. This suggests that the practice of copying code from other
contracts may introduce vulnerabilities and propagate the spread of vulnerable
contracts.

In future work, we would like to investigate other ways protocols can be used
to reduce gas usage. For instance, if two fields are used in different states,
it may be possible to use the same storage for both fields, eliminating the
cost of reserving an extra storage spot. We would also be interested in
studying protocol-related asset bugs, where a contract’s assets are lost or
can no longer be accessed due to a state transition. We would also like to
reduce the false positive rates on our StateOptimizationDetector and
AccessLockDetector to make them more useful to Solidity developers. We hope
this work can inform and motivative future research into the design of
languages and type systems for smart contracts.

## References

  * (1)
  * Beckman et al. (2011) Nels E. Beckman, Duri Kim, and Jonathan Aldrich. 2011\.  An Empirical Study of Object Protocols in the Wild. In _ECOOP 2011 – Object-Oriented Programming_ , Mira Mezini (Ed.). Springer Berlin Heidelberg, Berlin, Heidelberg, 2–26.
  * Buterin and Swende (2021) Vitalik Buterin and Martin Swende. 2021.  EIP-3529: Reduction in refunds.  https://eips.ethereum.org/EIPS/eip-3529
  * Chen et al. (2017) Ting Chen, Xiaoqi Li, Xiapu Luo, and Xiaosong Zhang. 2017\.  Under-optimized smart contracts devour your money. In _2017 IEEE 24th International Conference on Software Analysis, Evolution and Reengineering (SANER)_. IEEE, 442–446.
  * Coblenz et al. (2020) Michael Coblenz, Reed Oei, Tyler Etzel, Paulette Koronkevich, Miles Baker, Yannick Bloem, Brad A. Myers, Joshua Sunshine, and Jonathan Aldrich. 2020.  Obsidian: Typestate and Assets for Safer Blockchain Programming.  _ACM Trans. Program. Lang. Syst._ 42, 3, Article 14 (Nov. 2020), 82 pages.  https://doi.org/10.1145/3417516
  * Durieux et al. (2020) Thomas Durieux, João F. Ferreira, Rui Abreu, and Pedro Cruz. 2020.  Empirical Review of Automated Analysis Tools on 47,587 Ethereum Smart Contracts. In _Proceedings of the ACM/IEEE 42nd International Conference on Software Engineering_ (Seoul, South Korea) _(ICSE ’20)_. Association for Computing Machinery, New York, NY, USA, 530–541.  https://doi.org/10.1145/3377811.3380364
  * Elsden et al. (2018) Chris Elsden, Arthi Manohar, Jo Briggs, Mike Harding, Chris Speed, and John Vines. 2018\.  _Making Sense of Blockchain Applications: A Typology for HCI_.  Association for Computing Machinery, New York, NY, USA, 1–14.  https://doi.org/10.1145/3173574.3174032
  * Feist et al. (2019) Josselin Feist, Gustavo Grieco, and Alex Groce. 2019\.  Slither: A Static Analysis Framework for Smart Contracts.  _2019 IEEE/ACM 2nd International Workshop on Emerging Trends in Software Engineering for Blockchain (WETSEB)_ (May 2019).  https://doi.org/10.1109/wetseb.2019.00008
  * Foundation (2021) Ethereum Foundation. 2021\.  Common Patterns.  https://docs.soliditylang.org/en/v0.8.6/common-patterns.html#state-machine
  * IBM (2019) IBM. 2019.  Blockchain for Supply Chain - IBM Blockchain.  https://www.ibm.com/blockchain/industries/supply-chain
  * Kalra et al. (2018) Sukrit Kalra, Seep Goel, Mohan Dhawan, and Subodh Sharma. 2018\.  Zeus: Analyzing safety of smart contracts.. In _NDSS_. 1–12.
  * Kondo et al. (2020) Masanari Kondo, Gustavo Oliva, Zhen Jiang, Ahmed E. Hassan, and Osamu Mizuno. 2020.  Code Cloning in Smart Contracts: A Case Study on Verified Contracts from the Ethereum Blockchain Platform.  _Empirical Software Engineering_ 25 (11 2020).  https://doi.org/10.1007/s10664-020-09852-5
  * Luu et al. (2016) Loi Luu, Duc-Hiep Chu, Hrishi Olickel, Prateek Saxena, and Aquinas Hobor. 2016.  Making Smart Contracts Smarter. In _Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security_ (Vienna, Austria) _(CCS ’16)_. Association for Computing Machinery, New York, NY, USA, 254–269.  https://doi.org/10.1145/2976749.2978309
  * Moosavi and Clark (2021) Mahsa Moosavi and Jeremy Clark. 2021.  Lissy: Experimenting with on-chain order books.  arXiv:2101.06291 [cs.CR]
  * Ortner and Eskandari (2021) Martin Ortner and Shayan Eskandari. 2021.  Smart Contract Sanctuary.  (July 2021).  https://github.com/tintinweb/smart-contract-sanctuary
  * Popper (2016) Nathaniel Popper. 2016\.  A Hacking of More Than $50 Million Dashes Hopes in the World of Virtual Currency.  https://www.nytimes.com/2016/06/18/business/dealbook/hacker-may-have-removed-more-than-50-million-from-experimental-cybercurrency-project.html
  * Rozen (2021) Jacob Rozen. 2021\.  The ridiculously high cost of Gas on Ethereum.  https://coingeek.com/the-ridiculously-high-cost-of-gas-on-ethereum/
  * Schrans et al. (2018) Franklin Schrans, Susan Eisenbach, and Sophia Drossopoulou. 2018\.  Writing Safe Smart Contracts in Flint. In _Conference Companion of the 2nd International Conference on Art, Science, and Engineering of Programming_ (Nice, France) _(Programming’18 Companion)_. Association for Computing Machinery, New York, NY, USA, 218–219.  https://doi.org/10.1145/3191697.3213790
  * Strom and Yemini (1986) Robert E. Strom and Shaula Yemini. 1986.  Typestate: A programming language concept for enhancing software reliability.  _IEEE Transactions on Software Engineering_ SE-12, 1 (1986), 157–171.  https://doi.org/10.1109/TSE.1986.6312929
  * Tikhomirov et al. (2018) Sergei Tikhomirov, Ekaterina Voskresenskaya, Ivan Ivanitskiy, Ramil Takhaviev, Evgeny Marchenko, and Yaroslav Alexandrov. 2018\.  SmartCheck: Static Analysis of Ethereum Smart Contracts. In _2018 IEEE/ACM 1st International Workshop on Emerging Trends in Software Engineering for Blockchain (WETSEB)_. 9–16.
  * Wood et al. (2014) Gavin Wood et al. 2014\.  Ethereum: A secure decentralised generalised transaction ledger.  _Ethereum project yellow paper_ 151, 2014 (2014), 1–32.