# General theory of specific binding:
insights from a genetic-mechano-chemical protein model

John M. McBride<EMAIL_ADDRESS>Center for Soft and Living Matter,
Institute for Basic Science, Ulsan 44919, South Korea  Jean-Pierre Eckmann
<EMAIL_ADDRESS>Département de Physique Théorique and Section de
Mathématiques, University of Geneva, Geneva, Switzerland  Tsvi Tlusty
<EMAIL_ADDRESS>Center for Soft and Living Matter, Institute for Basic
Science, Ulsan 44919, South Korea  Departments of Physics and Chemistry, Ulsan
National Institute of Science and Technology, Ulsan 44919, South Korea

###### Abstract

Proteins need to selectively interact with specific targets among a multitude
of similar molecules in the cell. But despite a firm physical understanding of
binding interactions, we lack a general theory of how proteins evolve high
specificity. Here, we present such a model that combines chemistry, mechanics
and genetics, and explains how their interplay governs the evolution of
specific protein-ligand interactions. The model shows that there are many
routes to achieving molecular discrimination – by varying degrees of
flexibility and shape/chemistry complementarity – but the key ingredient is
precision. Harder discrimination tasks require more collective and precise
coaction of structure, forces and movements. Proteins can achieve this through
correlated mutations extending far from a binding site, which fine-tune the
localized interaction with the ligand. Thus, the solution of more complicated
tasks is enabled by increasing the protein size, and proteins become more
evolvable and robust when they are larger than the bare minimum required for
discrimination. The model makes testable, specific predictions about the role
of flexibility and shape mismatch in discrimination, and how evolution can
independently tune affinity and specificity. Thus, the proposed theory of
specific binding addresses the natural question of “why are proteins so big?”.
A possible answer is that molecular discrimination is often a hard task best
performed by adding more layers to the protein.

## Introduction

Proteins are the main molecular workforce inside cells, and the tasks they
perform invariably rely on specific, short-range interactions. The cell is
filled with thousands of molecular species, some differing by only a single
atom. Yet somehow, most proteins can specifically bind to only a few select
species. 1, 2, 3 Despite knowing much about pairwise binding, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14 we understand much less the many-body problem of how
proteins selectively bind to targets, while avoiding interactions with
similar, but non-cognate, molecules. 10 Such undesirable interactions can, at
best, lead to inefficiencies through inhibition, 15 and at worst, result in
aggregation 16, inaccurate translation, or cross-talk between signals. 17

Unwanted interactions can be minimized by designing mismatch between ligands
and binding pockets, 18, 19 such that the energetic cost of deformation allows
the protein to sift the target from similar, non-cognate ligands – a form of
‘conformational proofreading’. 20, 21 Recent work has suggested that residues
not directly involved in binding may play a role in discrimination via
allostery. 22, 23, 24, 25 Here, taking inspiration from the many experimental
examinations of molecular discrimination by proteins, such as enzymes, tRNA
synthetases, transcription factors, and antibodies , 26, 27, 1, 28, 29, 2, 30,
31, 32, 33 we propose a simple yet general theory of specific binding.

We study the evolution of discrimination by proteins using a genetic-mechano-
chemical model of binding. We find that, although the discrimination problem
is in general difficult, it has many possible solutions: shape mismatch,
chemical complementarity, and flexibility can all be manipulated in various
ways to tune interaction specificity. The important common thread is that it
requires precision and coordination – for example, just the right amount of
shape mismatch to fit a given flexibility. We show that residues distant from
the binding site enable this fine-tuning of mechanical deformation upon
binding – demonstrating that discrimination is the outcome of concerted,
collective interactions throughout the protein. Thus, larger proteins benefit
from having more degrees of freedom, allowing them to solve harder
discrimination challenges. Furthermore, larger proteins are more evolvable and
robust since their set of possible functional sequences is larger and more
connected. We further explain the mechanisms through which affinity and
specificity can be tuned independently, and discuss the role of flexibility
and entropy. Altogether, this simple model shines light on the difficult
problem of how proteins achieve such superlative molecular discrimination. At
the same time, by linking protein size to a property as ubiquitous as binding,
we offer a possible answer to the question, “why are proteins so big?”. 34, 35

## Results

Figure 1:  Mechano-chemical (MeCh) model of molecular recognition.  A: A
protein is modelled as a 2-d spring network with a chemical binding site
(three interaction sites, shown in pink); the $\Lambdaup$-shaped ligands have
three interaction sites, and are uniquely defined by an angle $\theta$ (three
examples are shown). B: Binding free energy gap in a rigid protein
$\upDelta\upDelta G_{{{\mathrm{L}}}{{\mathrm{S}}}}=\upDelta
G_{{\mathrm{L}}}-\upDelta G_{{\mathrm{S}}}$ as a function of the size of the
larger ligand $\theta_{{\mathrm{L}}}$, for various ligand size differences
$\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}$; minima are indicated with
dotted lines. C: $\upDelta\upDelta G_{{{\mathrm{L}}}{{\mathrm{S}}}}$ as a
function of the spring constant $K$ for
$\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}=$$$, for
$$$\leq\theta_{{\mathrm{L}}}\leq$$$. D: $\upDelta\upDelta
G_{{{\mathrm{L}}}{{\mathrm{S}}}}$ as a function of $\theta_{{\mathrm{L}}}$ for
$\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}=$$$, for $10\leq
K\leq${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}^{2}$$.

The models.  We introduce three protein models, starting with a bare-bones
model and then adding archetypal protein features, so that we achieve a
graduated understanding of how proteins evolve molecular discrimination. We
first introduce the basic mechano-chemical protein model (MeCh), a homogeneous
spring network with chemical binding sites that bind to a ligand (Fig. 1A).
This basic model shows how specific binding depends on an interplay of: shape
(mis)match between ligands and the protein binding pocket, protein
flexibility, chemical binding strength, and entropy. In the second, genetic
model (G-MeCh) we introduce a protein sequence which determines the stiffness
of individual springs, and the chemical binding strength of binding sites.
This allows us to study how specific binding is achieved in a model where
discrete changes occur via sequence mutations. In the final variant of the
model (G-MeCh-S), we allow the equilibrium structure of the protein to change
as a function of sequence, and we examine proteins of different sizes. This
model facilitates the study of protein evolution, and the effect of protein
size, via quantification of evolvability and mutational robustness.

Mechano-chemical (MeCh) model of molecular recognition.  The basic mechano-
chemical model is a 2-dimensional elastic network, 36, 37, 38 with
${n_{{\mathrm{A}}}}$ amino acids, $a_{i}$, arranged on a hexagonal lattice,
and ${n_{{\mathrm{B}}}}=3$ chemical binding sites, each with binding energy
constant $\varepsilon$ (Fig. 1A). In this coarse-grained representation, one
may envision the nodes as groups of tightly-connected amino acids that have
highly-correlated motion. 39 Likewise, the binding sites can be thought of as
a subdivision of a binding site into the three most salient units of amino
acids and functional groups. 40

The protein has a $\Lambdaup$-shaped binding pocket, described by an angle of
$\theta_{0}=$$$; opposite the protein is a set of $\Lambdaup$-shaped ligands
with three binding sites, which are uniquely described by an angle in the
range $$$\leq\theta\leq$$$. Shape mismatch is thus a function of only the
ligand shape, and varies along a single dimension,
$\upDelta\theta_{0}=\theta-\theta_{0}$. Regarding ligands, we use the
following terminology: Ligands that do not fit in the binding pocket
($\theta>$$$) are called ‘fat’, otherwise ‘thin’. When comparing two ligands,
the larger one is denoted ligand L ($\theta_{{\mathrm{L}}}$), and the smaller
one is denoted ligand S ($\theta_{{\mathrm{S}}}$). Similarly, we refer to the
cognate ligand ($\theta_{{\mathrm{C}}}$, i.e. the target), and the non-cognate
ligand(s) ($\theta_{{\mathrm{NC}}}$, the functionally undesirable one).

For each protein-ligand pair, we calculate the free energy of binding
$\upDelta G$ (i.e., the binding affinity) as the sum over three contributions:
the deformation energy ${\mathcal{E}}$, the chemical energy ${\mathcal{C}}$,
and change in entropy upon binding $\upDelta S$,

$\upDelta G={\mathcal{E}}+{\mathcal{C}}-T\upDelta S~{},$ (1)

where $T$ is temperature. Deformation energy is calculated as

${\mathcal{E}}=\frac{1}{2}\sum_{\langle i,j\rangle}K_{ij}(\Delta
r_{ij}-\ell_{ij})^{2}~{},$ (2)

where the summation $\langle i,j\rangle$ is over all ordered pairs of amino
acids $a_{i}$ and $a_{j}$ connected by a bond, $K_{ij}$ is the spring constant
of the bond, $\Delta r_{ij}$ is its length, and $\ell_{ij}$ its equilibrium
length. The initial configuration is not deformed, ${\mathcal{E}}=0$, since
all distances are equal to the equilibrium bond lengths, $\Delta
r_{ij}=\ell_{ij}$. In the MeCh model all bonds have identical lengths,
$\ell_{ij}=\ell=$1\text{\,}\mathrm{n}\mathrm{m}$$, and spring constants,
$K_{ij}=K$.

The chemical binding energy is given by

${\mathcal{C}}=\sum^{{n_{{\mathrm{B}}}}}_{i=1}-\varepsilon_{i}\,e^{-|{\mathbf{r}}^{b}_{i}-{\mathbf{r}}^{p}_{i}|^{2}/\sigma^{2}}~{},$
(3)

where $\varepsilon_{i}$ is the energy scale of binding locus $i$ of the
ligand, ${\mathbf{r}}^{p}_{i}$ and ${\mathbf{r}}^{b}_{i}$ are the positions,
respectively, of the amino acids $a_{j}$ at the binding pocket and the ligand
binding locus $i$, and we set the length scale of the interaction at
$\sigma=$0.3\text{\,}\mathrm{nm}\text{/}$$; interaction sites are uniquely
paired, so the energy is summed only over the ${n_{{\mathrm{B}}}}=3$ binding
pairs. In the MeCh model the energy scale is the same for all binding sites,
$\varepsilon_{i}=\varepsilon$.

The entropy change by binding $\upDelta S$ in Eq.(1) is the logarithm of the
relative change in the volume of the configuration space accessible by thermal
fluctuations of the protein. Thus, reduction in the magnitude of thermal
motion results in entropy loss. We decompose the binding entropy into two
terms, $\upDelta S=\upDelta S_{\mathrm{conf}}+\upDelta S_{0}$, where $\upDelta
S_{\mathrm{conf}}$ is the change in conformational entropy, 41, 42, 43 which
depends on the elastic network topology and spring stiffness, and $\upDelta
S_{0}$ accounts for other contributions to entropy that are not captured
directly by the model (e.g., release of frustrated solvent, ligand
conformational change, protein entropy change due to plastic deformation). 44,
45, 46, 47, 48, 49, 24

We calculate $\upDelta S_{\mathrm{conf}}$ for the elastic network by creating
stiff bonds of strength $K_{{\mathrm{\Lambdaup}}}$ between the protein and
ligand binding sites (see Methods). Standard normal mode analysis shows that
the resulting entropy change is the sum over the variation in the logarithms
of the mode energies $\lambda_{n}$ before and after binding, $\upDelta
S_{\mathrm{conf}}=-{\textstyle\frac{1}{2}}\sum_{n}{\Delta\ln\lambda_{n}}$. By
constraining the motion, stiffening typically increases the mode energies,
$\Delta\lambda_{n}\geq 0$, and binding therefore induces entropy loss,
$\upDelta S_{\mathrm{conf}}\leq 0$. For a homogeneous spring network
($K_{ij}=K$), this entropy change is well approximated as $\upDelta
S_{\mathrm{conf}}\approx-2\ln(K/K_{{\mathrm{\Lambdaup}}})+1.5$ (SI Fig. 1).

We compute the minimum binding free energy using a gradient descent algorithm
with some relevant physical constraints. Calculating this binding free energy
for multiple ligands gives us the binding free energy gap between any pair of
two ligands, $L$ and $S$, $\upDelta\upDelta
G_{{{\mathrm{L}}}{{\mathrm{S}}}}=\upDelta G_{{\mathrm{L}}}-\upDelta
G_{{\mathrm{S}}}$ (or between cognate and non-cognate ligands,
$\upDelta\upDelta G=\upDelta G_{{\mathrm{NC}}}-\upDelta G_{{\mathrm{C}}}$),
which we use as a measure of molecular discrimination and specificity.

Recognition via lock-and-key binding.  As a starting point, we examine the
limiting case of a completely rigid protein (i.e., $K=\infty$), which
corresponds to lock-and-key binding. 50 As no deformation occurs upon binding
and no conformational entropy is lost (if we set $\upDelta S_{0}=0$), the
binding free energy only depends on the mismatch
$\upDelta\theta_{0}=\theta-\theta_{0}$ of the shape of the wedge and the
ligand:

$\upDelta G_{\rm rigid}={\mathcal{C}}(r)=\begin{cases}-\varepsilon,&\text{if
}\theta>\theta_{0}\\\
-\varepsilon\left(2+e^{-r^{2}/\sigma^{2}}\right),&\text{if
}\theta\leq\theta_{0}\\\ \end{cases}~{},$

where
$r=2\ell\sin{{\textstyle\frac{1}{2}}\upDelta\theta_{0}}\approx\ell\upDelta\theta_{0}$
is the gap between the third interacting pair of the binding site and the
ligand. This is because fat ligands can only interact via one site, and thin
ligands minimize binding energy by fully binding to two sites, and partially
binding to the third.

In this limiting case, the best binding gap, $\upDelta\upDelta
G=2\varepsilon$, is achieved via steric exclusion: when one ligand is a
perfect match for the binding site ($\theta_{0}-\theta_{{\mathrm{S}}}=0$), and
the other ligand is fat, for any
$\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}=\theta_{{\mathrm{L}}}-\theta_{{\mathrm{S}}}$
(Fig. 1B). If both ligands are thin, the largest binding gap is much lower
($\upDelta\upDelta G=\varepsilon$), 33 and can only be obtained for
sufficiently dissimilar ligands
($\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}\geq$$$). If the two thin
ligands are similar, the optimal binding gap is only achievable with some
mismatch between the cognate ligand and the protein
($\theta_{{\mathrm{L}}}<\theta_{0}$). We thus see that binding with a rigid
protein is in principle a feasible strategy for molecular discrimination.
However, even in this hypothetical case, shape mismatch is often necessary to
promote binding of the cognate ligand while avoiding binding of the non-
cognate.

Figure 2:  Effect of energy and entropy on binding and discrimination in the
MeCh model.  A: Binding free energy for ligand L vs. ligand S;
$\varepsilon=$8\text{\,}\mathrm{k}\mathrm{T}$$, $-T\upDelta
S_{0}=$10\text{\,}\mathrm{k}\mathrm{T}$$,
$K_{{\mathrm{\Lambdaup}}}=${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}$^{2}$.
Colored regions indicate “loose discrimination” ($\upDelta G_{{\mathrm{C}}}<0$
and $\upDelta\upDelta G\geq$2\text{\,}\mathrm{k}\mathrm{T}$$), or “strict
discrimination” ($\upDelta G_{{\mathrm{C}}}<0$, $\upDelta G_{{\mathrm{NC}}}>0$
and $\upDelta\upDelta G\geq$2\text{\,}\mathrm{k}\mathrm{T}$$). Binding free
energy is shown for the proteins from Fig. 1C-D
($\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}=$$$). B: Binding free energy
gap as a function of spring constant $K$, and chemical energy constant
$\varepsilon$;
$K_{{\mathrm{\Lambdaup}}}=${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}$^{2}$,
$\theta_{{\mathrm{C}}}=$$$, $\theta_{{\mathrm{NC}}}=$$$. Optimal $K$ for a
given $\varepsilon$ is indicated by the orange dashed line. C: Regions that
satisfy the conditions for loose and strict discrimination for different
values of $-T\upDelta S_{0}$. D: Affinity ($\upDelta G$) versus specificity
($\upDelta\upDelta G$) for the parameter space sampled in (B-C).

Preferential binding by flexible proteins exploits ligand-protein shape
mismatch.  Rigidity in proteins is limited by the nature of the non-covalent
bonds holding them together. To explore the effect of rigidity on
discrimination, we first set wide bounds on the spring constants ($10\leq
K\leq${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}^{2}$$),
and later narrow these to more biologically-plausible values of spring
constants. 51, 52, 53 When $K$ is low, binding is not very specific as the
protein can perfectly match any ligand with little deformation energy (Fig.
1C). As $K$ increases, specificity improves, until the protein is so rigid
that neither cognate nor non-cognate ligands can induce the protein to deform.
The maximum binding free energy gap for fat ligands ($\upDelta\upDelta
G=0.84\varepsilon$, $\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}=$$$,
$K=${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}^{2}$$)
is considerably lower than what can be achieved by rigid lock-and-key binding
($\upDelta\upDelta G=2\varepsilon$).

For thin ligands, on the other hand, the maximum gap ($\upDelta\upDelta
G=0.49\varepsilon$, $\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}=$$$,
$K\sim$3\text{\cdot}{10}^{3}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}^{2}$$)
is double (Fig. 1C) of what can be achieved by a completely rigid protein
($\upDelta\upDelta G=0.24\varepsilon$,
$\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}=$$$). In general, higher
specificity is still achieved with greater rigidity across biologically-
relevant values of $K$; but this requires ever-greater shape match between the
protein and the cognate ligand, as $\upDelta\upDelta G$ becomes a steeper
function of ligand shape (Fig. 1D). This coupling between flexibility and
shape mismatch highlights the need for precise concerted control over both
protein structure and dynamics.

Discrimination is more difficult than recognition.  The binding energy gap is
a key determinant of molecular discrimination, but binding will only occur
spontaneously if the corresponding free energy change $\upDelta G$ is
negative, which depends on binding entropy. 41, 42 Taking this into account,
we now formally define molecular discrimination: ‘loose’ discrimination is
defined by $\upDelta\upDelta G\geq$2\text{\,}\mathrm{k}\mathrm{T}$$, and
$\upDelta G_{{\mathrm{C}}}<$0\text{\,}\mathrm{k}\mathrm{T}$$; we use
$2\text{\,}\mathrm{k}\mathrm{T}$ as a reasonable threshold that corresponds to
a sevenfold difference in binding affinity; ‘strict’ discrimination is defined
by further requiring $\upDelta
G_{{\mathrm{NC}}}>$0\text{\,}\mathrm{k}\mathrm{T}$$.

In Fig. 2A we replot the data from Fig. 1C-D in terms of discrimination,
finding that specific binding is substantially more difficult than recognition
alone. Even for our permissive threshold of $2\text{\,}\mathrm{k}\mathrm{T}$,
only a fraction of cases result in discrimination. The disparity between the
ability to distinguish between sets of fat ($\upDelta
G_{{\mathrm{S}}}<\upDelta G_{{\mathrm{L}}}$; upper left) versus thin
($\upDelta G_{{\mathrm{L}}}<\upDelta G_{{\mathrm{S}}}$; bottom right) ligands
illustrates the utility of steric exclusion as a discrimination mechanism
(discrimination is possible over a large range of $\upDelta S_{0}$ for fat
ligands, but if the ligands are thin then it is difficult to achieve
discrimination).

Affinity and specificity are correlated.  To understand the role of mechanical
and chemical binding energy in discrimination, we now vary the chemical energy
constant so that the chemical binding energy is within a reasonable biological
range ($1\leq\varepsilon\leq$25\text{\,}\mathrm{k}\mathrm{T}$$) . 54, 55, 56,
57, 58 Once again, increasing $K$ aids specificity, up until a point where the
protein is too rigid to accommodate either ligand (Fig. 2B). Beyond this
point, increases in $K$ must also be accompanied by higher $\varepsilon$, so
that there is sufficient chemical driving force for mechanical deformation.
$\varepsilon$ sets an upper bound on specificity ($\upDelta\upDelta G$), which
is achieved at some optimal $K$, and this optimal $K$ is proportional to
$\varepsilon$. As a result, affinity ($\upDelta G_{{\mathrm{C}}}$) and
specificity ($\upDelta\upDelta G$) are naturally correlated within the
parameter space (Fig. 2D, Pearson’s $r=0.38$, $p\ll 0.05$).

Entropy can be used to modulate affinity.  Binding entropy is a key
determinant of molecular discrimination (Fig. 2C). 55, 59, 60, 58, 43 For the
same discrimination task as in Fig. 2B, only loose discrimination is possible
for $-T\upDelta S_{0}=$0\text{\,}\mathrm{k}\mathrm{T}$$ since both ligands
bind to the protein; for this particular case, in the regions with a large
binding energy gap, $$3$<-T\upDelta
S_{\mathrm{conf}}<$10\text{\,}\mathrm{k}\mathrm{T}$$, is too low to avoid
binding to the non-cognate ligand. As we increase the entropy cost $-T\upDelta
S_{0}$, loose discrimination becomes more difficult since stronger chemical
bonds are needed for binding. However, increasing $-T\upDelta S_{0}$ also
enables strict discrimination, eventually to the point where the entropic cost
is too high for any ligand to bind. Since we have assumed that entropic cost
applies equally to cognate and non-cognate ligands, it does not affect
$\upDelta\upDelta G$; tuning entropy thus offers a way to decouple affinity
and specificity (Fig. 2D).

Figure 3:  Linking mechanics and binding chemistry to genetic sequence in the
G-MeCh model.  A: Amino acids are either $\mathbf{w}$ (blue, forming weak
bonds), or $\mathbf{s}$ (orange, forming strong bonds), as indicated by the
colour of the spheres; edge width indicates the three possible values of bond
strength of the bonds as determined by the sequence; chemical energy is
similarly determined by sequence as $\mathbf{w}_{c}$, or $\mathbf{s}_{c}$ as
indicated by the colour of the haloes at the interaction sites. Sequence is
read from left-to-right, from bottom-to-top (as indicated by the numbers).
B-C: Sequence correlation for all loose-discrimination solutions (B), and
solutions with $\upDelta\upDelta G>3.6$ kT (C). D-E: Mean sequence correlation
versus distance between residues; shaded region indicates the standard
deviation from the mean sequence correlation from $N_{S}$ random sequences.
For this example, the parameters are: $K_{\mathrm{s}}=${10}^{3}$$,
$K_{\mathrm{m}}=$400$$, $K_{\mathrm{w}}=$100$$,
$K_{{\mathrm{\Lambdaup}}}=${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}^{2}$$;
$\varepsilon_{\mathrm{s}}=6$,
$\varepsilon_{\mathrm{w}}=$3.88\text{\,}\mathrm{k}\mathrm{T}$$;
$\theta_{{\mathrm{S}}}=$$$, $\theta_{{\mathrm{L}}}=$$$; $T\upDelta
S_{0}=$0\text{\,}\mathrm{k}\mathrm{T}$$.

Sequence-dependent model (G-MeCh).  We have so far ignored the quintessential
feature of proteins – proteins are heteropolymers, composed of distinct amino
acids whose sequence is encoded in genes that are subject to evolution. We now
expand our model to include these two important components – heterogeneity in
flexibility, and change via discrete mutations – by coupling the mechanics and
chemistry to the protein sequence.

Thus, we examine a model protein consisting of ${n_{{\mathrm{A}}}}=13$ amino
acid letters, and ${n_{{\mathrm{B}}}}=3$ binding letters, using a 2-letter
alphabet for each (Fig. 3A). Amino acids are either $\mathbf{w}$ or
$\mathbf{s}$, and the elastic bond between each pair of neighbours,
$a_{i},a_{j}$, depends on their identities $K_{ij}(a_{i},a_{j})$: bonds are
either strong ($\mathbf{s}{-}\mathbf{s}$, with spring constant
$K_{\mathrm{s}}$), medium ($\mathbf{s}{-}\mathbf{w}$ or
$\mathbf{w}{-}\mathbf{s}$, $K_{\mathrm{m}}$), or weak
($\mathbf{w}{-}\mathbf{w}$, $K_{\mathrm{w}}$). Similarly, chemical binding
strength is determined by the letters $\mathbf{w}_{c}$ and $\mathbf{s}_{c}$,
respectively, resulting in weak ($\varepsilon_{\mathrm{w}}$) and strong bonds
($\varepsilon_{\mathrm{s}}$). The change in conformational entropy is
typically greater for more flexible proteins, and is incorporated in the model
by comparing the conformational entropy of the spring network before and after
binding (Methods).

Sequence variation reveals epistasis.  For a representative set of parameters,
we calculate the binding free energy $\upDelta G$ for all possible
$2^{{n_{{\mathrm{A}}}}+{n_{{\mathrm{B}}}}}=2^{16}=$65,536$$ sequences, and
identify those sequences that result in discrimination. On average, there is
no bias towards a particular amino acid at any position (SI Fig. 2), with the
exception of strong chemical binding at the tip of the $\Lambdaup$-shaped
binding pocket (position 15). However, detailed examination of sequence
correlation (Methods) reveals substantial, non-random patterns of epistasis
(Fig. 3B). There are positive correlations between the chemical binding sites
and the central amino acids in rows 2 and 3 – increasing the stiffness of the
intramolecular bonds made by these amino acids allows for stronger chemical
binding. One could anticipate this, given the correlations between
$\varepsilon$, $K$ and $\upDelta\upDelta G$ found in Fig. 2B. Thus, the subset
of solutions with $\upDelta\upDelta G>$3.6\text{\,}\mathrm{k}\mathrm{T}$$
(Fig. 3C) contains only those sequences that have taken advantage of this
positive correlation – all sequences have strong chemical bonds, and the
central amino acids typically have strong bonds (SI Fig. 2).

Far-away residues enable fine-tuning of protein binding.  All epistasis
visible in Fig. 3C is the outcome of negative correlations in the sequence.
For two positions that negatively covary, when one position goes from weak to
strong (or vice versa), the other position tends to do the opposite. In this
way protein flexibility is fine-tuned, to achieve a nearly-optimal amount of
deformation at the binding site. This is most apparent in the correlation
between all of the central seven amino acids. This can be explained by the
existence of many different sequences that can encode the same open-close
motion at the binding site (SI Fig. 3). Note, however, that this effective
fine-tuning of the protein mechanics is not limited to the binding site and is
evident throughout the full length of the protein, between pairs of residues
at distances up to $2.5\text{\,}\mathrm{nm}\text{/}$ (Fig. 3D-E).

Figure 4:  The G-MeCh-S model.  A: Protein models of different size. B:
Protein sequence is mapped to perturbations in the protein structure according
to an interaction matrix. If the interaction is positive, neighbouring amino
acids experience an attractive force with displacement magnitude $\delta
r=$0.1\text{\,}\mathrm{nm}\text{/}$$ and vice versa (black arrows); the
overall displacement is the average of the interactions with all neighbouring
amino acids. The open black circles are the original positions, and the
coloured circles show both the identities of the amino acids and the new
equilibrium configuration. C-D: Evolvability (C) and robustness (D) (both
normalized by sequence length) as a function of the angle of the cognate
ligand, $\theta_{{\mathrm{C}}}$, for easy, medium, and hard discrimination
tasks. Colours indicate protein models of different size. Green circles are
large to highlight that they overlap with purple circles. Symbols are only
shown for ligands for which solutions were found. Parameters are:
$K_{\mathrm{s}}=${10}^{3}$$, $K_{\mathrm{m}}=$400$$, $K_{\mathrm{w}}=$100$$,
$K_{{\mathrm{\Lambdaup}}}=${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}^{2}$$;
$\varepsilon_{\mathrm{s}}=6$,
$\varepsilon_{\mathrm{w}}=$3.88\text{\,}\mathrm{k}\mathrm{T}$$; $T\upDelta
S_{0}=$0\text{\,}\mathrm{k}\mathrm{T}$$.

Structural perturbation model (G-MeCh-S).  Residues far from the binding site
give the protein extra degrees of freedom that can be used to fine-tune
binding. We thus study how these degrees of freedom are related to the ability
to evolve functional sequences that are robust to mutations, by modelling
proteins of different size (Fig. 4A). To study protein evolution, we map
genotype to fitness via a simple measure $F$ of binding specificity,21

$F=\frac{e^{-\upDelta G_{{\mathrm{C}}}}-\sum_{i}{e^{-\upDelta
G_{{\mathrm{NC}}}^{i}}}}{1+e^{-\upDelta
G_{{\mathrm{C}}}}+\sum_{i}{e^{-\upDelta G_{{\mathrm{NC}}}^{i}}}}~{},$ (4)

where $\upDelta G_{{\mathrm{C}}}$ is the binding energy of the cognate ligand,
and the non-cognate ligands have binding energies $\upDelta
G_{{\mathrm{NC}}}^{i}$ (all energies are in kT units). The performance
measure, $-1\leq F\leq 1$, is the difference between the binding probability
of the cognate ligand and the sum of probabilities to bind any other
competitor, and sequences are considered fit if $F>0$.

The structure of real proteins varies with their sequence, generating
elaborately rugged fitness landscapes, with numerous maxima and basins of
attraction. 61, 62 While the G-MeCh model generates many fitness maxima, they
are rather shallow, forming large basins of nearly-optimal configurations.
This is because the shape of the G-MeCh protein is sequence-independent,
rendering the fitness landscape unrealistically smooth. Therefore, we
introduce the G-MeCh-S model, where the protein equilibrium structure is
allowed to change depending on sequence variation – a more realistic
approximation. The resulting affinity, specificity, and fitness landscapes
become much more rugged, containing multiple basins of attraction (SI Fig. 4).
As a protein’s shape varies, it evolves to switch preferences among competing
ligands. Note that we do not model protein folding, but instead account for
sequence variation resulting in small perturbations to the native structure
(Fig. 4B, Methods).

We study the effect of protein size on the evolution of proteins that can
discriminate a target ligand from many non-cognate ligands. We successively
increase task difficulty by increasing the number of non-cognate ligands. For
each task, the target ligand is picked from a set of ligands,
$\theta\in\\{$$,$$,\ldots,100^{\circ}\\}$. For each cognate ligand
$\theta_{{\mathrm{C}}}$ we designate its neighbours as competing non-cognate
ligands (i.e., those which satisfy
$|\theta_{{\mathrm{C}}}-\theta_{{\mathrm{NC}}}|\leq\phi$). This allows us to
consider discrimination tasks of variable difficulty: easy, $\phi=$$$, medium,
$\phi=$$$, and hard, $\phi=$$$, such that in the easiest task there are at
most 2 non-cognate ligands, and in the hardest task there are at most 6 non-
cognate ligands. In this way, when the target ligand has high shape mismatch,
there are fewer non-cognate ligands, so task difficulty is not strongly
dependent on $\theta_{{\mathrm{C}}}$. In principle, one could have also varied
task difficulty by changing $\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}$,
but the present approach has practical advantages (SI Fig. 4).

Larger proteins are more evolvable and robust.  We quantify evolvability as
the maximum number of mutations a fit protein can undergo while still
remaining fit ($F\geq 0$). Conversely, we quantify robustness as the minimum
number of mutations a fit protein can undergo before it becomes unfit. We find
that larger proteins are able to solve discrimination tasks for a greater
range of cognate ligands, and longer sequences are more evolvable and robust,
63, 64, 65 even after controlling for sequence length (Fig. 4C-D, Methods).
Increasing protein size gives diminishing returns which eventually saturate,
indicating that finite degrees of freedom are sufficient to achieve a
maximally evolvable and robust protein. For more difficult tasks, proteins of
different size exhibit a larger difference in evolvability and robustness, and
larger proteins are needed to reach the saturation point (Fig. 4C-D, SI Fig.
5). Thus, the required degrees of freedom depend on the difficulty of the
discrimination problem.

Figure 5:  Phase diagram of affinity and specificity for the MeCh model.  A-B:
Binding affinity (A) and specificity (B) as a function of protein flexibility,
and shape mismatch between cognate ligand and binding pocket
($\varepsilon=$8\text{\,}\mathrm{k}\mathrm{T}$$,
$|\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}|=$$$, $-T\upDelta S_{0}=0$,
$K_{{\mathrm{\Lambdaup}}}=${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}^{2}$$).
In the top plot (B) the cognate ligand is smaller than the non-cognate, and
the binding pocket is open to the ligand; in the bottom plot (B) the cognate
ligand is larger than the non-cognate, and the binding pocket is narrow
compared to the ligand. The optimal mismatch for a given flexibility is shown
by the black line. Cyan arrows (A) indicate proposed experimental trajectories
of antibodies during affinity maturation. 66 Purple arrows (B) indicate
proposed experimental trajectories of enzyme evolution (top – esterase, bottom
– cytochrome P450), and red arrows indicate predicted trajectories of
suggested experiments. 30, 67 C: Scaling of optimal mismatch $r^{*}$ with
chemical energy $\varepsilon$ and stiffness $K$. Fits (lines) to model results
(symbols) are shown for open (top) and closed (bottom) binding pockets, and
for $\varepsilon=2$ and $20$ kT. Symbols tend to overlap.

Tuning affinity and specificity via the interplay of energy, entropy, and
shape.  Proteins have an optimal affinity range for cognate ligand(s), 68 and
require some degree of specificity for functionality. The MeCh model explains
how these can be controlled by varying protein flexibility, shape, chemistry,
and entropy; this is illustrated in a phase diagram of affinity and
specificity (Fig. 5A-B). We study a smaller parameter space for the G-MeCh
model, and find qualitatively similar results (SI Fig. 6).

Affinity ($\upDelta G$) depends non-linearly on chemical bond strength, shape
mismatch, pocket geometry (open / closed), flexibility, and entropy (Fig. 5A).
The number, and strength of potential chemical bonds (represented in the model
by $\varepsilon$) sets an upper limit to affinity (if $-\upDelta S_{0}\geq 0$;
see “Limitations” section for counter-examples not covered here), which is
attained when a ligand fits perfectly into a rigid protein pocket. For rigid
proteins, shape mismatch sharply reduces affinity, but still results in zero-
deformation partial binding, which is stronger in proteins with open binding
pockets. Flexibility can both increase and decrease affinity: flexible
proteins have higher conformational entropy, which may result in higher
binding entropy cost and weaker affinity; 69 but if there is some shape
mismatch between the ligand and the binding pocket, flexibility increases
affinity by enabling the protein to adopt an optimal conformation for binding,
as clearly demonstrated by our model (Fig. 5A).

The specificity $\upDelta\upDelta G$ for a pair of cognate and non-cognate
ligands differing by $\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}=$$$ is
depicted in Fig. 5B. Maximum specificity for a given $\varepsilon$ is achieved
by matching flexibility and shape mismatch (Fig. 5B), so that deformation
energy is maximized (SI Fig. 7). 40 Proteins need to be just flexible enough
for the cognate ligand to bind, but not the non-cognate ligand. Too much
flexibility (reduces deformation energy) and both will bind; 70 too rigid and
both will partially bind with negligible deformation. The optimal deformation
for a given flexibility is found to scale linearly with the chemical binding
energy $\varepsilon$ (Fig. 5C), such that the overall shape of the specificity
phase diagram is not affected by changes to $\varepsilon$ (SI Fig. 8). We will
later demonstrate how these findings can be recapitulated by a simplified
phenomenological model.

The model predicts how entropy may affect affinity and specificity: The
primary effect of entropy in the MeCh-model is to modulate affinity via
changes in conformational entropy $\upDelta S_{\mathrm{conf}}$ as the
flexibility varies (Fig. 5A). This stiffening effect will be similar for
competing ligands and thereby will impact specificity only weakly (although
for sufficiently dissimilar ligands, this contribution of entropy to
specificity may be more significant). Thus, by selectively tuning entropy, one
can modulate affinity almost independently of specificity (Fig. 2D). In
practice, this might be difficult; for example, flexibility simultaneously
affects binding energy and conformational entropy. 59 There is more potential
for independent control via the other components to binding entropy, $\upDelta
S_{0}$. Examples of this include the use of intrinsically-disordered regions:
through coupled folding and binding; 71, 72 or through a disordered region far
from the binding site, which is affected allosterically. 49 Other routes to
decoupling specificity and affinity include control over oligomeric
complexation, 73 or interactions with solvent. 46, 47

Phenomenological model.  The general nature of the specificity can be
illustrated by a simple phenomenological model (Methods): the free energy is
taken as a function of a single coordinate $r$, which represents the gap
between a pair of interacting loci of the binding site and the ligand. The
initial gap $r_{0}$, before the protein deforms, is proportional to the
mismatch, $r_{0}\approx\ell\upDelta\theta_{0}$. The elastic energy is a
function of the deformation $u=r_{0}-r$, while the chemical energy is a
function of the gap $r$, to which we add binding entropy to obtain the free
energy (Eq.(1)), $\upDelta
G(r,r_{0})={\mathcal{E}}(u)+{\mathcal{C}}(r)-T\upDelta S~{}$. This protein-
ligand system reaches an equilibrium configuration, $r=\bar{r}$ when the
elastic and chemical forces counterbalance each other
$-{\mathcal{E}}^{\prime}(u)+{\mathcal{C}}^{\prime}(r)=0$ (where ′ is the
derivative).

The specificity $\upDelta\upDelta G$ in discriminating between a cognate
ligand and a competitor, whose initial gaps differ by $\delta r_{0}$, is the
change in $\upDelta G$ along the equilibrium line, $\bar{r}(r_{0})$. Analysis
shows that the specificity is proportional to the forces,
${\mathcal{D}}=\upDelta\upDelta G/\delta
r_{0}\approx{\mathcal{E}}^{\prime}(u)={\mathcal{C}}^{\prime}(r)$ (where
${\mathcal{D}}$ is the “specific specificity”, or specificity per unit of
difference between ligands). Therefore, maximal specificity is obtained
exactly for the mismatch $r_{0}^{\ast}$ and gap $r^{\ast}$ at which the
opposing elastic and chemical forces are both maximal, at the inflection
point, ${\mathcal{C}}^{\prime\prime}(r)=0$, where the system is most sensitive
to small shape differences (see details in Methods). This simple result
explains why the line of optimal specificity follows the line of maximal
deformation energy in our model (SI Fig. 7).

We apply this general theory to MeCh model, and consider the standard Hookean
elastic energy ${\mathcal{E}}(u)={\textstyle\frac{1}{2}}K_{\mathrm{eff}}u^{2}$
(Eq.(2)). The whole protein is expected to be softer than a single bond, with
an effective spring constant $K_{\mathrm{eff}}=\alpha K$ ($\alpha\leq 1$). As
for the chemical energy (Eq.(3)), we assume that two loci of the
$\Lambdaup$-shaped ligand are in contact with the (open) binding pocket and
contribute $-2\varepsilon$ to the free energy, which together with the energy
of the third locus gives
${\mathcal{C}}(r)=-\varepsilon(2+e^{-r^{2}/\sigma^{2}})$. The inflection point
of this chemical energy is located at $r^{\ast}=\sigma/\sqrt{2}$, such that
the optimal mismatch is

$\frac{r_{0}^{\ast}}{\sigma}=2^{-{\textstyle\frac{1}{2}}}+(2/e)^{\textstyle\frac{1}{2}}\left(\frac{\varepsilon}{K_{\mathrm{eff}}\sigma^{2}}\right)\approx
0.71+\frac{0.85}{\alpha}\left(\frac{\varepsilon}{K\sigma^{2}}\right)~{},$ (5)

with a corresponding optimal “specific specificity”,
${\mathcal{D}}^{\ast}=\upDelta\upDelta G^{\ast}/\delta r_{0}\approx
0.85(\varepsilon/\sigma)$.

We fit Eq.(5) to the MeCh model predictions for the optimal mismatch (Fig.
5C). For the case of an open binding pocket ($\upDelta\theta_{0}>0$), we find
an intercept of $0.5-0.6$ and effective softening $0.57<\alpha<0.71$, in broad
agreement with the phenomenological model. For a narrow pocket, we find
qualitatively similar behaviour: due to steric exclusion, less mismatch is
needed for a maximally rigid protein (i.e., the intercept is closer to zero);
the slope of the $r_{0}^{\ast}/\sigma$ line is approximately double in this
case, since the reduction in affinity due to steric exclusion is roughly
twofold ($2\varepsilon$) compared to the open pocket case (Fig. 1B). The
phenomenological model (Eq.(5)) elucidates two findings of the detailed
protein models: (i) Mismatch is needed even for extremely rigid proteins (Fig.
5B). (ii) The optimal mismatch depends on a trade-off between chemical binding
energy and mechanical energy, as demonstrated in Fig. 2B where optimal
specificity is achieved at a constant ratio $\varepsilon/K$. In principle,
this simple fit will also apply to the G-MeCh model, but there
$K_{\mathrm{eff}}(\\{K_{ij}\\})$ has to be extracted for each protein (SI Fig.
6).

## Discussion

Through analysis of a conceptually simple, yet multi-faceted model of
molecular discrimination, we start to understand the mechanisms, and
evolution, of molecular discrimination.

How difficult is molecular discrimination?  It depends on the context. For
example, proteins that bind to nucleic acids or lipids must find their target
out of a dizzying gallery of lookalikes. 74, 28, 75 Enzymes need to be able to
release their product after catalysis, 76 which can be more (e.g.,
isomerization reactions) or less (e.g., proteolysis) challenging depending on
how similar the product is to the substrate. In principle, the difficulty,
$\Psi$, can be expressed as a function of the set of ligands,
$\\{\Lambdaup\\}$, and the required degree of specificity,

$\Psi\left[\Lambdaup^{{\mathrm{C}}},\upDelta
G_{{\mathrm{C}}},\Lambdaup^{{\mathrm{NC}}}_{1},\upDelta\upDelta
G_{1},\ldots~{},\Lambdaup^{{\mathrm{NC}}}_{i},\upDelta\upDelta
G_{i}\right]~{},$

where $\upDelta G_{{\mathrm{C}}}$ is the optimum binding affinity of the
cognate ligand, and $\upDelta\upDelta G_{i}$ is the required minimal binding
energy gap between the cognate ligand $\Lambdaup^{{\mathrm{C}}}$ and non-
cognate ligand $i$, $\Lambdaup^{{\mathrm{NC}}}_{i}$; the degree of specificity
required for each non-cognate ligand depends on the in vivo concentration, 70
and the cost of incorrect binding.

The difficulty $\Psi$ is a non-linear function, which we expect to generally
increase with: the required specificity ($\upDelta\upDelta G_{i}$), the number
of similar, non-cognate ligands $\Lambdaup^{{\mathrm{NC}}}_{i}$, and how easy
it is to distinguish them from the cognate ligand $\Lambdaup^{{\mathrm{C}}}$.
It is easier to distinguish: if the cognate ligand is smaller than the non-
cognate, as it allows steric exclusion of the larger ligand (Fig. 1B); if the
cognate ligand can form more energetic bonds (e.g., an extra hydroxyl group
effectively increases $\varepsilon$); 77 or has distinct chemical differences
(e.g., positive versus negative charge) to the non-cognate ligands. 78 We lack
a general, robust method of quantifying $\Psi$; creating a metric of ligand
discriminability (not the same as similarity 79, 80) would be immensely useful
for, e.g., predicting sites for specific inhibition of proteins. For now, we
propose that an ad hoc approximation for $\Psi$ could be taken to be the
degrees of freedom needed to achieve functional discrimination (Fig. 6A).

Figure 6:  Theory of biomolecular discrimination.  A: More difficult
discrimination tasks require more degrees of freedom, which in proteins
corresponds to longer sequences. Sequences that are longer than the minimum
necessary size are more robust and evolvable. B: Larger proteins have more
degrees of freedom with which to fine-tune structure and dynamics at the
binding site, since mutations at distal residues can have small yet
significant effects.

Molecular discrimination by a hypothetical mechano-chemical machine.  Lock-
and-key binding is most specific, but only under very strict conditions (Fig.
1B): The non-cognate ligands must be larger than the cognate ligand (steric
exclusion), the cognate ligand must perfectly match the binding site, and both
protein and ligand must be rigid. Otherwise, some mismatch between cognate
ligand and binding site is always needed to optimise specificity (Fig. 1D) –
such ‘conformational proofreading’ (i.e., an energetic penalty) minimises
binding to the non-cognate ligand(s). 20, 21, 81, 82 The path to better
discrimination still lies with higher rigidity, but increasing control over
shape mismatch is needed, since deviations from the optimal mismatch are
tolerated less and less as rigidity increases (Fig. 1D). Ultimately, the key
feature of good discrimination is precision: the right amount of flexibility,
coupled with the right shape, results in optimal deformation.

Molecular discrimination by proteins.  Proteins, as _genetic_ -mechano-
chemical machines, have some inherent features that strongly constrain their
molecular discrimination ability. 77, 83 Most interactions between amino acids
are non-covalent, so proteins cannot be very rigid. Moreover, protein geometry
is limited to the topology of a folding chain, composed of discrete units
(amino acids) of approximately $1\text{\,}\mathrm{nm}\text{/}$ in size, so
perfect shape match is practically impossible. Finally, evolution advances in
discrete steps, mutation, deletions and insertions – not via continuous
tuning.

In light of these constraints, mutations very close to a binding site are
bound to have a large effect on the flexibility and geometry of the binding
site. 84 On the other hand, mutations further from the binding site can have
ever smaller effects, 85, 86 enabling fine-tuning of mechanics (Fig. 3), 87,
88, 89, 90, 91, 92 and structure. 93, 94 For example, in a separate recent
study of ours, we studied the effect of single mutations on structure (using
proteins in the protein data bank, and AlphaFold), finding that structural
perturbations are felt up to $2\text{\,}\mathrm{n}\mathrm{m}$ away from the
mutated residue. 95 This is consistent with our assertion that far-away
mutations can influence molecular discrimination.

These observations, taken together with the finding that good discrimination
necessitates precision, lead us to propose that larger proteins – as they have
more degrees of freedom (potential mutations) – can achieve better
discrimination through finer control over protein dynamics and structure at
the binding site (Fig. 6B). For a given discrimination task there is a minimum
protein size (Fig. 4C-D), but such an efficient protein may be difficult to
evolve. As protein size grows beyond the bare minimum, there are ever more
sequences capable of solving the problem, which results in sequences that are
evolvable and robust (Fig. 4C-D). The proposed theory thus predicts that
proteins have a lower bound on the size required to achieve discrimination,
and that they will be larger than this in order to be evolvable (Fig. 6A). It
is not clear how close proteins get to this minimum, but it would be more
likely in prokaryotes since they have greater efficiency requirements.

Examining and predicting experimental trends.  Affinity and specificity vary
non-monotonically with $\upDelta\theta_{{{\mathrm{L}}}{{\mathrm{S}}}}$, $K$,
and $\varepsilon$, so there are no general, unidirectional trends. Still, we
can discuss trends found in the experimental literature in the context of Fig.
5A-B, and offer explanations that are consistent with our model and lead to
testable predictions.

Many studies report that germline antibodies are flexible, and become more
rigid in a process known as affinity maturation. 96, 97 A recent in silico
directed evolution study corroborated this, and also showed that some
antibodies first become more flexible, before later becoming more rigid. 66
According to our model (cyan arrows, Fig. 5A), the former case should occur
when antibodies have close to optimal shape mismatch (low
$|\upDelta\theta_{0}|$); the latter case should occur when antibodies have
high shape mismatch (high $|\upDelta\theta_{0}|$). This prediction can be
tested by measuring shape (and chemistry; note we operationalized mismatch
through shape, but shape and chemistry are inseparable in real molecules)
mismatch: one can use measures based on static structure, 98, 99, 100 but for
better results one should calculate how often the optimal binding
configuration is sampled in the free antibody, using molecular dynamics
simulations. 101, 102, 103, 104

Our model supports the notion that ancestral enzymes were both flexible and
promiscuous, 105 as flexibility is often anti-correlated with specificity in
our model. However, we find that the correlation between flexibility and
specificity depends on shape and chemical binding energy (Fig. 5B). This can
be illustrated by two case studies that relate conflicting accounts of the
role of flexibility in enzyme promiscuity: Flexibility and promiscuity are
correlated in a group of 57 human cytochrome P450 (CYT) enzymes, 106 while the
opposite trend is observed in a group of 147 esterases. 67 In the case of CYT,
we know that the binding pocket is quite small ($\upDelta\theta_{0}<0$), so we
can infer that the proteins fall to the right of the optimal line (bottom
purple arrow, Fig. 5B). Thus, our model predicts that as rigidity increases,
both specificity and affinity will eventually decrease as the protein will be
too stiff to deform (bottom red arrow, Fig. 5B). This can be tested by
increasing CYT stiffness via directed evolution. 107 In contrast to CYT, the
esterases have open active sites ($\upDelta\theta_{0}>0$), and we know that
promiscuity is correlated with the volume of the active site, 78, 30 so we can
infer that the proteins fall to the left of the optimal line (top purple
arrow, Fig. 5B). In this case, further increasing esterase flexibility via
directed evolution should reveal that there is an optimal range of flexibility
where specificity is maximized (top red arrow, Fig. 5B).

Our model explains how affinity and specificity can be either positively 108,
109 or negatively 110 correlated within a set of proteins, depending on how
they differ in shape mismatch and flexibility (SI Fig. 7). For example,
affinity and specificity are positively correlated in two cases: when proteins
differ along optimal line (black line, Fig. 5A-B), or orthogonal to the
optimal line. In the former case, deformation energy decreases when affinity
and specificity increase; in the latter case, deformation energy increases
when affinity and specificity increase (SI Fig. 7). This prediction may be
tested by studying the transcription factor Pho4, where increased binding
affinity of the cognate CACGTG nucleotide sequence was found to improve
discrimination of the cognate over the non-cognate CACGTT sequence in 210
variants. 94 These variants can be studied using molecular dynamics
simulations, where deformation energy in our model is analogous to the change
in internal energy of a protein upon binding 111, 112, 113.

A challenge in testing these predictions is the vast amount of data that is
needed, since one needs to measure multiple dimensions for a combination of
proteins and ligands. However, multiple methods can characterise shape
mismatch, flexibility, chemical bond energy, deformation energy, and entropy.
We advocate combining molecular dynamics simulations (which can characterise
flexibility and calculate deformation energy) with high-throughput experiments
(which can measure binding kinetics), and to develop methods to control for
orthogonal effects such as differences in protein stability or foldability. 93
Existing public data sets from previous experiments present a facile
opportunity in this regard. 1, 29, 30, 2, 106

Molecular discrimination by aminoacyl-tRNA synthetases.  To evaluate the
theory that protein size depends on task difficulty, we need to know both the
relevant non-cognate ligands and necessary binding specificity. In the case of
aminoacyl-tRNA synthetases (ARSs), we know the relevant ligands (the 20-30
proteogenic and non-proteogenic amino acids present in cells), and that they
have similar in vivo concentrations and similar costs associated with missense
mutations. This presents a natural control, such that task difficulty is
reduced to a question of discriminability between cognate and non-cognate
ligands. This is still difficult to evaluate, but we can start by using the
available experimental data on pairwise binding specificity of ARSs. We can
rationalise that one out of a pair is easier to recognize if it is smaller by
a methyl group (steric exclusion), or has an extra hydroxyl group (can form
more high energy bonds). Thus, it is difficult to discriminate: threonine from
serine (minus one methyl), isoleucine from valine (minus one methyl),
phenylalanine from tyrosine (plus one hydroxyl), and alanine from serine (plus
one hydroxyl). 33

Comparing ARSs of these pairs, we find that the ARS of the easier-to-recognize
ligand has greater specificity and, with the exception of Val-Ile, they are
also smaller (Fig. 7). This exception may be due to the difficulty in general
in distinguishing between many aliphatic amino acids, a point illustrated by
the fact that these ARSs all have post-transfer editing domains. 114
Furthermore, when we compare ARS enzymes to non-ARS enzymes that also act on
amino acids (but with lower specificity requirements), the non-ARS enzymes
tend to be considerably smaller than ARSs (SI Fig. 9). 33 These findings
support the theory that protein size is related to the difficulty of the
discrimination task. We expect that ARSs are an exemplary class with which to
further evaluate theories on protein specificity.

Figure 7:  Pairwise selectivity – ratio of cognate vs non-cognate $K_{\rm
M}/k_{\rm cat}$ – of aminoacyl-tRNA synthases (ARSs), and ARS complex size, of
pairs of similar amino acid ligands, where in each pair the labels are
coloured according to whether it should be easier (black) or more difficult
(red) to discern. Ligands differ by one methyl group (left, ligand with
additional $\mathrm{-CH_{3}}$ in red) or one hydroxyl group (right, additional
$\mathrm{-OH}$ in black). Differences are highlighted in diagrams (top) with
dotted circles.

Limitations and extensions.  By focusing on binding energy gap we have posed
the molecular discrimination as a thermodynamic problem. It is not immediately
clear whether a kinetics-focused approach would lead to the same conclusions.
For example, we find that increasing rigidity can increase specificity, but
from a kinetics point of view, one might expect the opposite. There is
evidence to suggest that flexibility aids formation of initial encounter
complexes, which reduces the rate of futile encounters. 9, 115, 13
Understanding the role of kinetics in molecular discrimination warrants a
separate, focused study.

We treat chemical binding strength and mismatch (encoded in shape) as
separable components in our model. This is useful for studying their relative
contributions to molecular discrimination, but molecules cannot be
deconstructed in this way. It may be possible, however, to describe the
differences in interactions between a protein and several ligands using some
type of principal components or reduced dimensions. In this way, one may
relate the model results to discrimination of real molecules. In a similar
vein, we have simplified our model by assuming that ligands are rigid. We
speculate that if ligands were allowed to move, this would produce results
consistent with a rigid ligand and a more flexibile protein. This is
consistent with our finding (via the phenomenological model) that many spring
constants can be reduced to one effective spring constant. Additionally, if
one ligand is more flexible than another, it may be captured by an ‘effective
mismatch’ term that takes into account a ligand’s ability to change shape. In
general, protein function is typically governed by a small number of effective
parameters, as manifested by drastic dimensional reduction in genotype-to-
phenotype map. 116 Our simple model demonstrates that specificity is a
function of the effective flexibility, shape mismatch, and local chemical
interaction, but further work on more detailed models is needed to understand
how many effective dimensions are needed to describe specificity.

Our model treats deformations as elastic, so it may not generalize well to
proteins that undergo plastic deformations upon binding. In fact, the
restructuring of intramolecular bonds can result in a gain in entropy44, 24,
or an increase in internal enthalpy. 117 Furthermore, a simplifying assumption
in elastic network models is that bonds are at their equilibrium lengths in
the ensemble-average protein structure. However, proteins are non-crystalline
matter where many bonds are “pre-stressed" due to geometrical frustration (the
inability to achieve equilibrium length in all bonds simultaneously). 118,
119. The frustration gives rise to a rugged energy landscape with myriad of
local minima. The model can be easily extended to incorporate frustration by
modeling a network with heterogeneous equilibrium lengths $\ell_{ij}$, and
considering a thermal ensemble of configurations. Accounting for the role of
plasticity and intrinsic disorder in biomolecular discrimination will require
a more complex model.

Why are proteins so big?  Proteins are large macromolecules. Large proteins
require large genomes, and thus slow transcription and translation.
Prokaryotic proteins, which prioritize small genomes and fast replication,
should benefit immensely from smaller proteins. We do see that proteins are on
average smaller in prokaryotes (312 residues) than in eukaryotes (441
residues), 120 but they are still quite large. Thus, there is some
indispensable functional reason for protein size. We have here proposed that
the difficulty of evolving proteins that can discriminate necessitates
proteins of a certain size. We now discuss some alternative factors affecting
protein size.

Stability is typically a base requisite for a functional protein, 35 but most
proteins are marginally stable – i.e., stability is not maximised, but rather
an acceptable level is reached. 83 Longer proteins can more easily fulfill
these requirements, 121 but short sequences are often sufficiently stable. 122
Catalytic activity in enzymes is extremely important, but comparable activity
levels can be found in much smaller organic catalysts. 123, 124 Large surface
area may be needed to make multiple interaction sites, whether for a single
molecule, 125 or multiple molecules. 34, 35, 96 While we propose that size is
determined by the difficulty of molecular discrimination, all of the above
constraints may also be relevant factors.

Data availability.  All code (models, analysis, figures), and some data is
available at github.com/jomimc/RecProt. Additional data can be shared upon
reasonable request.

Acknowledgements.  We thank Jacques Rougemont for discussions. This work was
supported by the Institute for Basic Science, Project Code IBS-R020-D1. JPE is
partially supported by the Fonds National Suisse Swissmap.

METHODS

Calculating binding entropy.  Binding will typically involve loss of entropy
due to constraining the internal fluctuations of the protein. To calculate
this change in conformational entropy upon binding, one starts by analyzing
the elastic energy of a spring network. It is given by
${\mathcal{E}}={\textstyle\frac{1}{2}}\bra{u}{\mathbf{H}}\ket{u}$ , where
$|u\rangle$ is the displacement vector (of size $d\cdot{n_{{\mathrm{A}}}}$)
and ${\mathbf{H}}$ is the elasticity matrix (for details see 25, 23). Spectral
decomposition of ${\mathbf{H}}$ gives
${\mathbf{H}}=\sum_{n}\lambda_{n}\outerproduct{u_{n}}{u_{n}}$, where
$\ket{u_{n}}$ are the normal modes and $\lambda_{n}$ are the eigenvalues.
There are ${\textstyle\frac{1}{2}}d(d+1)=3$ zero energy modes of translation
and rotation without deformation that we ignore. We can then express the
energy in normal mode coordinates, $\xi_{n}=\bra{u_{n}}\ket{u}$, as
${\mathcal{E}}={\textstyle\frac{1}{2}}\sum_{n}\lambda_{n}\xi_{n}^{2}$.

To estimate the entropy, we consider the partition function of the elastic
deformations,
${\mathcal{Z}}=\int\differential{\ket{u}}e^{-\beta{\mathcal{E}}}$, where
$\beta=1/(kT)$ is the inverse temperature. Expressed in normal coordinates, we
get
${\mathcal{Z}}=\prod_{n}{[\int\differential{\xi_{n}}\exp(-\beta\lambda_{n}\xi_{n}^{2}/2)]}=\prod_{n}[2\pi/(\beta\lambda_{n})]^{1/2}$,
where the product is taken over the modes for which $\lambda_{n}\neq 0$. The
entropy $S$ is obtained from the relation
$S=\ln{{\mathcal{Z}}}-\partialderivative*{\ln{{\mathcal{Z}}}}{\ln{\beta}}$, so
that $S={\textstyle\frac{1}{2}}\sum_{n}{[\ln(2\pi e)-\ln(\beta\lambda_{n})]}$.
When the protein binds, the modes and their eigenvalues will change. The
resulting entropy change is the sum over the changes in the logarithms of the
elastic mode energies before and after binding,

$\upDelta S=S_{\rm bound}-S_{\rm
unbound}=-{\textstyle\frac{1}{2}}\sum_{n}{\Delta\ln\lambda_{n}}~{}.$

Typically, the mode energies $\lambda_{n}$ increase or stay unchanged upon
binding since the motion is more constrained, and therefore the entropy is
reduced, $\upDelta S<0$, as the protein-ligand complex stiffens. 40 This
effect is calculated as follows: We add a new bond at the bottom of the
$\Lambdaup$ (opening of the binding site), and make the three bonds between
the amino acids of the binding site very rigid by increasing them to
$K_{{\mathrm{\Lambdaup}}}=$
${10}^{4}\text{\,}\mathrm{k}\mathrm{T}\mathrm{/}\mathrm{n}\mathrm{m}$2.

In the model of entropy used in this work, entropy decreases upon binding, and
flexible proteins lose more conformational entropy than rigid proteins. This
results from our choice of binding stiffness, assuming that binding stiffens
the pocket, $K_{{\mathrm{\Lambdaup}}}=1000\geq K$; and that
$K_{{\mathrm{\Lambdaup}}}$ is practically constant for all proteins. In
reality, flexible proteins will of course have more entropy, and thus more to
lose. Proteins are most often found to lose entropy upon binding, but many
proteins instead gain entropy due to allosteric conformational change. 44
Moreover, there are other contributions to entropy, such as solvent entropy
and ligand entropy, 126, 45, 48 that are not treated in our model (e.g.,
contributions of translational and rotational entropy do not depend on the
internal degrees of freedom and can be excluded).

Phenomenological model.  Consider a free energy $\upDelta G$ which is a
function of $r$, the remaining gap in the binding site (see text), and the
initial gap $r_{0}\approx\ell\upDelta\theta_{0}$. The elastic energy is a
function of the deformation $u=r_{0}-r$, while the chemical energy is a
function of the gap $r$, such that the overall free energy is $\upDelta
G(r,r_{0})={\mathcal{E}}(u)+{\mathcal{C}}(r)-T\upDelta S$, with the entropy
change $\upDelta S$. For a given mismatch $r_{0}$ the protein equilibrates at
a gap $r=\bar{r}$ when the elastic and chemical forces are equal in magnitude
and opposite, $(\partialderivative*{\upDelta
G}{r})_{r_{0}}=-{\mathcal{E}}^{\prime}(u)+{\mathcal{C}}^{\prime}(r)=0$. This
condition defines an equilibrium line in the $(r,r_{0})$-plane,
$M(\bar{r},r_{0})=0$. Taking the differential of the line,
$\variation{M}=-{\mathcal{E}}^{\prime\prime}(u)\variation{r_{0}}+({\mathcal{E}}^{\prime\prime}(u)+{\mathcal{C}}^{\prime\prime}(r))\variation{r}=0$,
we obtain its slope,
$\derivative*{r_{0}}{\bar{r}}=1+{\mathcal{C}}^{\prime\prime}(r)/{\mathcal{E}}^{\prime\prime}(u)$.
The binding free energy is the equilibrium value along the line
$\upDelta\bar{G}(r_{0})=\upDelta G(\bar{r},r_{0})$.

The specificity is the free energy difference between a cognate ligand and a
competitor, $\upDelta\upDelta G=\upDelta\bar{G}(r_{0}+\delta
r_{0})-\upDelta\bar{G}(r0)\approx{\mathcal{D}}\cdot\delta r_{0}$, where the
“specific specificity" (or “discriminabilty") is the specificity per shape
difference, ${\mathcal{D}}\equiv\derivative*{\upDelta\bar{G}}{r_{0}}$. Hence,
the specific specificity is exactly the elastic force,
${\mathcal{D}}=(\partialderivative*{\upDelta
G}{r_{0}})_{r}={\mathcal{E}}^{\prime}(u)$ (since along the equilibrium,
$(\partialderivative*{\upDelta G}{r})_{r_{0}}=0$). Maximal specificity is
therefore obtained for the mismatch $r_{0}^{\ast}$ at which the force is
maximal,$\derivative*{{\mathcal{E}}^{\prime}(u)}{r_{0}}={\mathcal{E}}^{\prime\prime}(u)(\derivative*{u}{r_{0}})={\mathcal{E}}^{\prime\prime}(u)(1-\derivative*{\bar{r}}{r_{0}})=0$.
It follows that maximal specificity is achieved when the slope of the
equilibrium curve is $\derivative*{r_{0}}{\bar{r}}=1$. From the equation of
the equilibrium line, we find that this optimum is the inflection point of the
chemical energy, ${\mathcal{C}}^{\prime\prime}(r)=0$, where the chemical
attractive force is also maximal. For chemical energies that lack an
inflection point, one searches for a global maximum of ${\mathcal{C}}(r)$.
There are also cases with multiple equilibrium points and discontinuities,
which need to be considered separately. In the text, we apply this general
result to the MeCh-model, and obtain the manifold of optimal $r_{0}^{\ast}$
mismatch as a function of model parameters $K$, and $\varepsilon$.

Calculating sequence covariance and correlation.  We follow the standard
procedure: One takes a set of ${n_{{\mathrm{S}}}}$ sequences $\mathbf{S}_{i}$
of length ${n_{{\mathrm{A}}}}+{n_{{\mathrm{B}}}}$, represented as vectors of
ones ($\mathbf{s}$) and zeros ($\mathbf{w}$), and subtracts from each the
average sequence vector $\mathbf{S}_{i}\to\mathbf{S}_{i}-\bar{\mathbf{S}}$
(all entries of $\mathbf{S}_{i}$ are between 0 and 1). The resulting sequences
are the rows of an
${n_{{\mathrm{S}}}}\times({n_{{\mathrm{A}}}}+{n_{{\mathrm{B}}}})$ matrix,
${\mathcal{M}}$. The covariance matrix is then
${\mathcal{Q}}=({n_{{\mathrm{S}}}}-1)^{-1}{\mathcal{M}}^{\rm T}{\mathcal{M}}$.
Finally, the matrix ${\mathcal{P}}$ of the Pearson correlation coefficients
between pairs of positions in the sequence (the “correlation matrix") is
${\mathcal{P}}=\mathrm{dg}{\mathcal{Q}}^{-1/2}{\mathcal{Q}}\,\mathrm{dg}{\mathcal{Q}}^{-1/2}$,
where $\mathrm{dg}{\mathcal{Q}}$ is the matrix of the diagonal elements of
${\mathcal{Q}}$ (i.e., the variances).

Connecting genes to structure.  We allow neighbouring $\mathbf{w}$ and
$\mathbf{s}$ amino acids to interact such that they either attract, repel, or
neither: attraction (repulsion) between two amino acids results in their
equilibrium positions moving closer together (further apart) by
$\delta{\mathbf{r}}=\alpha$ ($\delta{\mathbf{r}}=-\alpha$). We generate
interaction tables $\alpha(a_{i},a_{j})$ such that $\mathbf{w}{-}\mathbf{w}$,
$\mathbf{w}{-}\mathbf{s}$ / $\mathbf{s}{-}\mathbf{w}$, and
$\mathbf{s}{-}\mathbf{s}$ bonds, respectively result in one of the three
possible interactions. Thus, we can generate up to $15$ unique interaction
tables (reduced from a total of $3^{4}=81$ possible tables, by accounting for
symmetries). We show results for one set in Fig. 4; we verified that the
results do not depend on a particular set, and show results for another set in
SI Fig. 10. The new equilibrium position ${\mathbf{r}}_{i}$ of each amino acid
is then determined by the original equilibrium position
$\tilde{{\mathbf{r}}}_{i}$ plus the average over the displacements induced by
its $N_{j}$ neighbours,

${\mathbf{r}}_{i}=\tilde{{\mathbf{r}}}_{i}+\frac{1}{N_{j}}\sum_{j}^{N_{j}}\alpha(a_{i},a_{j})\,\hat{{\mathbf{r}}}_{ij}~{},$

where $j$ is the index of a bonded neighbour and
$\hat{{\mathbf{r}}}_{ij}=(\tilde{{\mathbf{r}}}_{j}-\tilde{{\mathbf{r}}}_{i})/|\tilde{{\mathbf{r}}}_{j}-\tilde{{\mathbf{r}}}_{i}|$
is a unit vector pointing from $\tilde{{\mathbf{r}}}_{i}$ to
$\tilde{{\mathbf{r}}}_{j}$. The equilibrium lengths of the springs are then
set to $\ell_{ij}=|{\mathbf{r}}_{i}-{\mathbf{r}}_{j}|$. An alternative
approach is to model mutational deformations as the response to linear
perturbations of the network’s Hamiltonian, 127, 128 as there is mounting
evidence of correlations between protein mechanics and structural evolution.
129, 23

The optimization algorithm.  For fixed coupling constants, we find a
configuration which optimizes energy. After moving the protein close to the
ligand, we use the L-BFGS method 130 to find a minimum. To model steric
repulsion, we restrict the positions of the amino acids so that they are
outside the sector defined by the ligand position and angle $\theta$. Two
post-processing checks precede the analysis: We require that the orientation
of each triangle in the protein maintains its orientation, i.e., that the
surface is not flipped, and the minimized energy must be at least as low as
that obtained by a completely rigid protein. We discard any results that do
not pass these tests.

Evolvability and robustness.  Evolvability is the ability of a population of
organisms to evolve new phenotypes. We measure this as the maximum number of
mutations a fit sequence can accumulate without having fitness less than zero.
Robustness is the ability of a sequence to mutate while retaining its original
function. We measure robustness as the minimum number of mutations needed for
a fit sequence to become non-fit. In practice, we calculated the edit distance
between fit sequences, and clustered them using a density-based clustering
algorithm (DBSCAN, implemented in the sklearn python module; $\mathrm{eps}=1$,
min_samples $=1$). 131 Evolvability of a sequence is then the maximum distance
within a cluster, while robustness is the minimum distance within a cluster.
For both measures we averaged over all fit sequences (according to Eq.(4)),
and normalized by dividing by the genome size.

## References

  * Nam  _et al._ (2012) H. Nam, N. E. Lewis, J. A. Lerman, D.-H. Lee, R. L. Chang, D. Kim, and B. O. Palsson, Network context and selection in the evolution to enzyme specificity, Science 337, 1101 (2012).
  * Piazza  _et al._ (2018) I. Piazza, K. Kochanowski, V. Cappelletti, T. Fuhrer, E. Noor, U. Sauer, and P. Picotti, A map of protein-metabolite interactions reveals principles of chemical communication, Cell 172, 358 (2018).
  * Copley (2020) S. D. Copley, The physical basis and practical consequences of biological promiscuity, Phys. Biol. 17, 051001 (2020).
  * Miller and Dill (1997) D. W. Miller and K. A. Dill, Ligand binding to proteins: The binding landscape model, Protein Sci. 6, 2166 (1997).
  * McCammon (1998) J. A. McCammon, Theory of biomolecular recognition, Curr. Opin. Struc. Biol. 8, 245 (1998).
  * Ma  _et al._ (1999) B. Ma, S. Kumar, C. J. Tsai, and R. Nussinov, Folding funnels and binding mechanisms, Protein Engineering, Design and Selection 12, 713 (1999).
  * Gohlke and Klebe (2002) H. Gohlke and G. Klebe, Approaches to the description and prediction of the binding affinity of small-molecule ligands to macromolecular receptors, Angewandte Chemie International Edition 41, 2644 (2002).
  * Mobley and Dill (2009) D. L. Mobley and K. A. Dill, Binding of small-molecule ligands to proteins: “what you see” is not always “what you get”, Structure 17, 489 (2009).
  * Wang  _et al._ (2013a) Q. Wang, P. Zhang, L. Hoffman, S. Tripathi, D. Homouz, Y. Liu, M. N. Waxham, and M. S. Cheung, Protein recognition and selection through conformational and mutually induced fit, P. Natl. Acad. Sci. Usa. 110, 20545 (2013a).
  * Kastritis and Bonvin (2013) P. L. Kastritis and A. M. J. J. Bonvin, On the binding affinity of macromolecular interactions: Daring to ask why proteins interact, J. Roy. Soc. . Interface 10, 20120835 (2013).
  * Weikl and Paul (2014) T. R. Weikl and F. Paul, Conformational selection in protein binding and function, Protein Sci. 23, 1508 (2014).
  * Stank  _et al._ (2016) A. Stank, D. B. Kokh, J. C. Fuller, and R. C. Wade, Protein binding pocket dynamics, Accounts Chem. Res. 49, 809 (2016).
  * Plattner  _et al._ (2017) N. Plattner, S. Doerr, G. De Fabritiis, and F. Noé, Complete protein–protein association kinetics in atomic detail revealed by molecular dynamics simulations and markov modelling, Nat. Chem. 9, 1005 (2017).
  * Wan  _et al._ (2020) S. Wan, A. P. Bhati, S. J. Zasada, and P. V. Coveney, Rapid, accurate, precise and reproducible ligand–protein binding free energy prediction, Interface Focus 10, 20200007 (2020).
  * Goldstein (1944) A. Goldstein, The mechanism of enzyme-inhibitor-substrate reactions: Illustrated by the cholinesterase-physostigmine-acetylcholine system, The Journal of general physiology 27, 529 (1944).
  * Mahler  _et al._ (2009) H.-C. Mahler, W. Friess, U. Grauschopf, and S. Kiese, Protein aggregation: Pathways, induction factors and analysis, J. Pharm. Sci-us. 98, 2909 (2009).
  * Olson and Marais (2000) M. F. Olson and R. Marais, Ras protein signalling, Semin. Immunol. 12, 63 (2000).
  * Rohs  _et al._ (2010) R. Rohs, X. Jin, S. M. West, R. Joshi, B. Honig, and R. S. Mann, Origins of specificity in protein-dna recognition, Annu. Rev. Biochem. 79, 233 (2010).
  * Sergeeva  _et al._ (2020) A. P. Sergeeva, P. S. Katsamba, F. Cosmanescu, J. J. Brewer, G. Ahlsen, S. Mannepalli, L. Shapiro, and B. Honig, Dip/dpr interactions and the evolutionary design of specificity in protein families, Nat. Commun. 11, 2125 (2020).
  * Savir and Tlusty (2007) Y. Savir and T. Tlusty, Conformational proofreading: The impact of conformational changes on the specificity of molecular recognition, Plos One 2, 1 (2007).
  * Savir and Tlusty (2013) Y. Savir and T. Tlusty, The ribosome as an optimal decoder: A lesson in molecular recognition, Cell 153, 471 (2013).
  * Rivoire (2019) O. Rivoire, Parsimonious evolutionary scenario for the origin of allostery and coevolution patterns in proteins, Phys. Rev. E 100, 032411 (2019).
  * Eckmann  _et al._ (2019) J. P. Eckmann, J. Rougemont, and T. Tlusty, Colloquium: Proteins: The physics of amorphous evolving matter, Rev. Mod. Phys. 91, 031001 (2019).
  * Wankowicz  _et al._ (2022) S. A. Wankowicz, S. H. de Oliveira, D. W. Hogan, H. van den Bedem, and J. S. Fraser, Ligand binding remodels protein side-chain conformational heterogeneity, eLife 11, e74114 (2022).
  * Dutta  _et al._ (2018) S. Dutta, J. P. Eckmann, A. Libchaber, and T. Tlusty, Green function of correlated genes in a minimal mechanical model of protein evolution, P. Natl. Acad. Sci. Usa. 115, E4559 (2018).
  * Sundberg and Mariuzza (2002) E. J. Sundberg and R. A. Mariuzza, Molecular recognition in antibody-antigen complexes, in  _Protein Modules and Protein-Protein Interaction_ , Vol. 61 (Academic Press, 2002) pp. 119–160.
  * Elias  _et al._ (2012) M. Elias, A. Wellner, K. Goldin-Azulay, E. Chabriere, J. A. Vorholt, T. J. Erb, and D. S. Tawfik, The molecular basis of phosphate discrimination in arsenate-rich environments, Nature 491, 134 (2012).
  * Jankowsky and Harris (2015) E. Jankowsky and M. E. Harris, Specificity and nonspecificity in rna–protein interactions, Nat. Rev. Mol. Cell Bio. 16, 533 (2015).
  * Huang  _et al._ (2015) H. Huang, C. Pandya, C. Liu, N. F. Al-Obaidi, M. Wang, L. Zheng, S. Toews Keating, M. Aono, J. D. Love, B. Evans, R. D. Seidel, B. S. Hillerich, S. J. Garforth, S. C. Almo, P. S. Mariano, D. Dunaway-Mariano, K. N. Allen, and J. D. Farelli, Panoramic view of a superfamily of phosphatases through substrate profiling, P. Natl. Acad. Sci. Usa. 112, E1974 (2015).
  * Martínez-Martínez  _et al._ (2018) M. Martínez-Martínez, C. Coscolín, G. Santiago, J. Chow, P. J. Stogios, R. Bargiela, C. Gertler, J. Navarro-Fernández, A. Bollinger, S. Thies, C. Méndez-García, A. Popovic, G. Brown, T. N. Chernikova, A. García-Moyano, G. E. K. Bjerga, P. Pérez-García, T. Hai, M. V. Del Pozo, R. Stokke, I. H. Steen, H. Cui, X. Xu, B. P. Nocek, M. Alcaide, M. Distaso, V. Mesa, A. I. Peláez, J. Sánchez, P. C. F. Buchholz, J. Pleiss, A. Fernández-Guerra, F. O. Glöckner, O. V. Golyshina, M. M. Yakimov, A. Savchenko, K. E. Jaeger, A. F. Yakunin, W. R. Streit, P. N. Golyshin, V. Guallar, M. Ferrer, and T. I. Consortium, Determinants and prediction of esterase substrate promiscuity patterns, Acs Chem. Biol. 13, 225 (2018).
  * Lite  _et al._ (2020) T. L. V. Lite, R. A. Grant, I. Nocedal, M. L. Littlehale, M. S. Guo, and M. T. Laub, Uncovering the basis of protein-protein interaction specificity with a combinatorially complete library, eLife 9, e60924 (2020).
  * Biswas and Thattai (2020) A. Biswas and M. Thattai, Promiscuity and specificity of eukaryotic glycosyltransferases, Biochem. Soc. T. 48, 891 (2020).
  * Tawfik and Gruic-Sovulj (2020) D. S. Tawfik and I. Gruic-Sovulj, How evolution shapes enzyme selectivity - lessons from aminoacyl-trna synthetases and other amino acid utilizing enzymes, The FEBS Journal 287, 1284 (2020).
  * Payens (1983) T. A. J. Payens, Why are enzymes so large?, Trends Biochem. Sci. 8, 46 (1983).
  * Srere (1984) P. A. Srere, Why are enzymes so big?, Trends Biochem. Sci. 9, 387 (1984).
  * Tirion (1996) M. M. Tirion, Large amplitude elastic motions in proteins from a single-parameter, atomic analysis, Phys. Rev. Lett. 77, 1905 (1996).
  * Chennubhotla  _et al._ (2005) C. Chennubhotla, A. J. Rader, L. W. Yang, and I. Bahar, Elastic network models for understanding biomolecular machinery: From enzymes to supramolecular assemblies, Phys. Biol. 2, S173 (2005).
  * López-Blanco and Chacón (2016) J. R. López-Blanco and P. Chacón, New generation of elastic network models, Curr. Opin. Struc. Biol. 37, 46 (2016).
  * Halabi  _et al._ (2009) N. Halabi, O. Rivoire, S. Leibler, and R. Ranganathan, Protein sectors: Evolutionary units of three-dimensional structure, Cell 138, 774 (2009).
  * Richard (2019) J. P. Richard, Protein flexibility and stiffness enable efficient enzymatic catalysis, J. Am. Chem. Soc. 141, 3320 (2019).
  * Frederick  _et al._ (2007) K. K. Frederick, M. S. Marlow, K. G. Valentine, and A. J. Wand, Conformational entropy in molecular recognition by proteins, Nature 448, 325 (2007).
  * Tzeng and Kalodimos (2012) S. R. Tzeng and C. G. Kalodimos, Protein activity regulation by conformational entropy, Nature 488, 236 (2012).
  * Sun  _et al._ (2017) Z. Sun, Y. N. Yan, M. Yang, and J. Z. H. Zhang, Interaction entropy for protein-protein binding, The Journal of Chemical Physics 146, 124124 (2017).
  * Fenley  _et al._ (2012) A. T. Fenley, H. S. Muddana, and M. K. Gilson, Entropy–enthalpy transduction caused by conformational shifts can obscure the forces driving protein–ligand binding, P. Natl. Acad. Sci. Usa. 109, 20006 (2012).
  * Dragan  _et al._ (2017) A. I. Dragan, C. M. Read, and C. Crane-Robinson, Enthalpy–entropy compensation: The role of solvation, Eur. Biophys. J. 46, 301 (2017).
  * Quiocho  _et al._ (1989) F. A. Quiocho, D. K. Wilson, and N. K. Vyas, Substrate specificity and affinity of a protein modulated by bound water molecules, Nature 340, 404 (1989).
  * Michel  _et al._ (2009) J. Michel, J. Tirado-Rives, and W. L. Jorgensen, Energetics of displacing water molecules from protein binding sites: Consequences for ligand optimization, J. Am. Chem. Soc. 131, 15403 (2009).
  * Peccati and Jiménez-Osés (2021) F. Peccati and G. Jiménez-Osés, Enthalpy–entropy compensation in biomolecular recognition: A computational perspective, ACS Omega 6, 11122 (2021).
  * Keul  _et al._ (2018) N. D. Keul, K. Oruganty, E. T. Schaper Bergman, N. R. Beattie, W. E. McDonald, R. Kadirvelraj, M. L. Gross, R. S. Phillips, S. C. Harvey, and Z. A. Wood, The entropic force generated by intrinsically disordered segments tunes protein function, Nature 563, 584 (2018).
  * Fischer (1894) E. Fischer, Einfluss der configuration auf die wirkung der enzyme, Ber. Dtsch. Chem. Ges. 27, 2985 (1894).
  * Atilgan  _et al._ (2001) A. R. Atilgan, S. R. Durell, R. L. Jernigan, M. C. Demirel, O. Keskin, and I. Bahar, Anisotropy of fluctuation dynamics of proteins with an elastic network model, Biophys. J. 80, 505 (2001).
  * Hinsen (2007) K. Hinsen, Structural flexibility in proteins: Impact of the crystal environment, Method. Biochem. Anal. 24, 521 (2007).
  * Riccardi  _et al._ (2009) D. Riccardi, Q. Cui, and G. N. Phillips, Application of elastic network models to proteins in the crystalline state, Biophys. J. 96, 464 (2009).
  * Schapira  _et al._ (1999) M. Schapira, M. Totrov, and R. Abagyan, Prediction of the binding energy for small molecules, peptides and proteins, J. Mol. Recognit. 12, 177 (1999).
  * Barratt  _et al._ (2005) E. Barratt, R. J. Bingham, D. J. Warner, C. A. Laughton, S. E. V. Phillips, and S. W. Homans, Van der waals interactions dominate ligand-protein association in a protein binding site occluded from solvent water, J. Am. Chem. Soc. 127, 11827 (2005).
  * Gallicchio  _et al._ (2010) E. Gallicchio, M. Lapelosa, and R. M. Levy, Binding energy distribution analysis method (bedam) for estimation of protein-ligand binding affinities, J. Chem. Theory Comput. 6, 2961 (2010).
  * Moal  _et al._ (2011) I. H. Moal, R. Agius, and P. A. Bates, Protein-protein binding affinity prediction on a diverse set of structures, Method. Biochem. Anal. 27, 3002 (2011).
  * Huang and Liu (2013) Y. Huang and Z. Liu, Do intrinsically disordered proteins possess high specificity in protein-protein interactions?, Chemistry - A European Journal 19, 4462 (2013).
  * Grünberg  _et al._ (2006) R. Grünberg, M. Nilges, and J. Leckner, Flexibility and conformational entropy in protein-protein binding, Structure 14, 683 (2006).
  * Chang  _et al._ (2007) C. A. Chang, W. Chen, and M. K. Gilson, Ligand configurational entropy and protein binding, P. Natl. Acad. Sci. Usa. 104, 1534 (2007).
  * Hietpas  _et al._ (2011) R. T. Hietpas, J. D. Jensen, and D. N. A. Bolon, Experimental illumination of a fitness landscape, P. Natl. Acad. Sci. Usa. 108, 7896 (2011).
  * de Visser and Krug (2014) J. A. G. M. de Visser and J. Krug, Empirical fitness landscapes and the predictability of evolution, Nat. Rev. Genet. 15, 480 (2014).
  * Levin  _et al._ (2009) K. B. Levin, O. Dym, S. Albeck, S. Magdassi, A. H. Keeble, C. Kleanthous, and D. S. Tawfik, Following evolutionary paths to protein-protein interactions with high affinity and selectivity, Nature Structural & Molecular Biology 16, 1049 (2009).
  * Tawfik and S. (2010) O. K. Tawfik and D. S., Enzyme promiscuity: A mechanistic and evolutionary perspective, Annu. Rev. Biochem. 79, 471 (2010).
  * Wagner (2013) A. Wagner,  _Robustness and Evolvability in Living Systems_ (Princeton university press, 2013).
  * Ovchinnikov  _et al._ (2018) V. Ovchinnikov, J. E. Louveau, J. P. Barton, M. Karplus, and A. K. Chakraborty, Role of framework mutations and antibody flexibility in the evolution of broadly neutralizing antibodies, eLife 7, e33038 (2018).
  * Nutschel  _et al._ (2021) C. Nutschel, C. Coscolín, B. David, D. Mulnaes, M. Ferrer, K. E. Jaeger, and H. Gohlke, Promiscuous esterases counterintuitively are less flexible than specific ones, J. Chem. Inf. Model. 61, 2383 (2021).
  * Popp  _et al._ (2020) A. P. Popp, J. Hettich, and J. C. M. Gebhardt, Transcription factor residence time dominates over concentration in transcription activation, bioRxiv 10.1101/2020.11.26.400069 (2020).
  * Cao  _et al._ (2016) H. Cao, Y. Huang, and Z. Liu, Interplay between binding affinity and kinetics in protein-protein interactions, Proteins. 84, 920 (2016).
  * Gade  _et al._ (2021) M. Gade, L. L. Tan, A. M. Damry, M. Sandhu, J. S. Brock, A. Delaney, A. Villar-Briones, C. J. Jackson, and P. Laurino, Substrate dynamics contribute to enzymatic specificity in human and bacterial methionine adenosyltransferases, JACS Au 1, 2349 (2021).
  * Wang  _et al._ (2013b) Y. Wang, X. Chu, S. Longhi, P. Roche, W. Han, E. Wang, and J. Wang, Multiscaled exploration of coupled folding and binding of an intrinsically disordered molecular recognition element in measles virus nucleoprotein, P. Natl. Acad. Sci. Usa. 110, E3743 (2013b).
  * Yang  _et al._ (2019) J. Yang, M. Gao, J. Xiong, Z. Su, and Y. Huang, Features of molecular recognition of intrinsically disordered proteins via coupled folding and binding, Protein Sci. 28, 1952 (2019).
  * Scheidt  _et al._ (2021) T. Scheidt, J. A. Carozza, C. C. Kolbe, F. A. Aprile, O. Tkachenko, M. M. J. Bellaiche, G. Meisl, Q. A. E. Peter, T. W. Herling, S. Ness, M. Castellana-Cruz, J. L. P. Benesch, M. Vendruscolo, C. M. Dobson, P. Arosio, and T. P. J. Knowles, The binding of the small heat-shock protein $\alpha$b-crystallin to fibrils of $\alpha$-synuclein is driven by entropic forces, P. Natl. Acad. Sci. Usa. 118, 10.1073/pnas.2108790118 (2021).
  * Kribelbauer  _et al._ (2019) J. F. Kribelbauer, C. Rastogi, H. J. Bussemaker, and R. S. Mann, Low-affinity binding sites and the transcription factor specificity paradox in eukaryotes, Annu. Rev. Cell Dev. Bi. 35, 357 (2019).
  * Contreras  _et al._ (2012) F. X. Contreras, A. M. Ernst, P. Haberkant, P. Björkholm, E. Lindahl, B. Gönen, C. Tischer, A. Elofsson, G. von Heijne, C. Thiele, R. Pepperkok, F. Wieland, and B. Brügger, Molecular recognition of a single sphingolipid species by a protein’s transmembrane domain, Nature 481, 525 (2012).
  * Rivoire (2020) O. Rivoire, Geometry and flexibility of optimal catalysts in a minimal elastic model, The Journal of Physical Chemistry B 124, 807 (2020).
  * Peracchi (2018) A. Peracchi, The limits of enzyme specificity and the evolution of metabolism, Trends Biochem. Sci. 43, 984 (2018).
  * Barrozo  _et al._ (2015) A. Barrozo, F. Duarte, P. Bauer, A. T. P. Carvalho, and S. C. L. Kamerlin, Cooperative electrostatic interactions drive functional evolution in the alkaline phosphatase superfamily, J. Am. Chem. Soc. 137, 9061 (2015).
  * Stephenson and Freeland (2013) J. D. Stephenson and S. J. Freeland, Unearthing the root of amino acid similarity, J. Mol. Evol. 77, 159 (2013).
  * Rácz  _et al._ (2018) A. Rácz, D. Bajusz, and K. Héberger, Life beyond the tanimoto coefficient: Similarity measures for interaction fingerprints, Journal of Cheminformatics 10, 48 (2018).
  * Savir and Tlusty (2010) Y. Savir and T. Tlusty, Reca-mediated homology search as a nearly optimal signal detection system, Mol. Cell 40, 388 (2010).
  * De Vlaminck  _et al._ (2012) I. De Vlaminck, M. van Loenhout, L. Zweifel, J. den Blanken, K. Hooning, S. Hage, J. Kerssemakers, and C. Dekker, Mechanism of homology recognition in dna recombination from dual-molecule experiments, Mol. Cell 46, 616 (2012).
  * Sikosek and Chan (2014) T. Sikosek and H. S. Chan, Biophysics of protein evolution and evolutionary protein biophysics, J. Roy. Soc. . Interface 11, 20140419 (2014).
  * Lopes  _et al._ (2010) A. Lopes, M. S. A. Busch, and T. Simonson, Computational design of protein-ligand binding: Modifying the specificity of asparaginyl-trna synthetase, J. Comput. Chem. 31, 1273 (2010).
  * Adams  _et al._ (2019) R. M. Adams, J. B. Kinney, A. M. Walczak, and T. Mora, Epistasis in a fitness landscape defined by antibody-antigen binding free energy, Cell Systems 8, 86 (2019).
  * Wang  _et al._ (2020) M. Wang, Z. Cang, and G. W. Wei, A topology-based network tree for the prediction of protein–protein binding affinity changes following mutation, Nature Machine Intelligence 2, 116 (2020).
  * Koenig  _et al._ (2017) P. Koenig, C. V. Lee, B. T. Walters, V. Janakiraman, J. Stinson, T. W. Patapoff, and G. Fuh, Mutational landscape of antibody variable domains reveals a switch modulating the interdomain conformational dynamics and antigen binding, P. Natl. Acad. Sci. Usa. 114, E486 (2017).
  * Otten  _et al._ (2018) R. Otten, L. Liu, L. R. Kenner, M. W. Clarkson, D. Mavor, D. S. Tawfik, D. Kern, and J. S. Fraser, Rescue of conformational dynamics in enzyme catalysis by directed evolution, Nat. Commun. 9, 1314 (2018).
  * Modi  _et al._ (2021) T. Modi, V. A. Risso, S. Martinez-Rodriguez, J. A. Gavira, M. D. Mebrat, W. D. Van Horn, J. M. Sanchez-Ruiz, and S. Banu Ozkan, Hinge-shift mechanism as a protein design principle for the evolution of $\beta$-lactamases from substrate promiscuity to specificity, Nat. Commun. 12, 1852 (2021).
  * Ose  _et al._ (2022) N. J. Ose, B. M. Butler, A. Kumar, I. C. Kazan, M. Sanderford, S. Kumar, and S. B. Ozkan, Dynamic coupling of residues within proteins as a mechanistic foundation of many enigmatic pathogenic missense variants, Plos Comput. Biol. 18, 1 (2022).
  * Ray  _et al._ (2021) D. Ray, L. Le, and I. Andricioaei, Distant residues modulate conformational opening in sars-cov-2 spike protein, P. Natl. Acad. Sci. Usa. 118, 10.1073/pnas.2100943118 (2021).
  * Karamitros  _et al._ (2022) C. S. Karamitros, K. Murray, B. Winemiller, C. Lamb, E. M. Stone, S. D’Arcy, K. A. Johnson, and G. Georgiou, Leveraging intrinsic flexibility to engineer enhanced enzyme catalytic activity, P. Natl. Acad. Sci. Usa. 119, e2118979119 (2022).
  * Markin  _et al._ (2021) C. J. Markin, D. A. Mokhtari, F. Sunden, M. J. Appel, E. Akiva, S. A. Longwell, C. Sabatti, D. Herschlag, and P. M. Fordyce, Revealing enzyme functional architecture via high-throughput microfluidic enzyme kinetics, Science 373, eabf8761 (2021).
  * Aditham  _et al._ (2021) A. K. Aditham, C. J. Markin, D. A. Mokhtari, N. DelRosso, and P. M. Fordyce, High-throughput affinity measurements of transcription factor and dna mutations reveal affinity and specificity determinants, Cell Systems 12, 112 (2021).
  * McBride  _et al._ (2022) J. M. McBride, K. Polev, V. Reinharz, B. A. Grzybowski, and T. Tlusty, Alphafold2 can predict structural and phenotypic effects of single mutations, bioRxiv 10.1101/2022.04.14.488301 (2022).
  * Van Regenmortel (2014) M. H. V. Van Regenmortel, Specificity, polyspecificity, and heterospecificity of antibody-antigen recognition, J. Mol. Recognit. 27, 627 (2014).
  * Mishra and Mariuzza (2018) A. K. Mishra and R. A. Mariuzza, Insights into the structural basis of antibody affinity maturation from next-generation sequencing, Frontiers in Immunology 9, 117 (2018).
  * Chen and Weng (2003) R. Chen and Z. Weng, A novel shape complementarity scoring function for protein-protein docking, Proteins. 51, 397 (2003).
  * Chen and Honig (2010) B. Y. Chen and B. Honig, Vasp: A volumetric analysis of surface properties yields insights into protein-ligand binding specificity, Plos Comput. Biol. 6, 1 (2010).
  * Yan and Huang (2019) Y. Yan and S. Y. Huang, Pushing the accuracy limit of shape complementarity for protein-protein docking, Bmc Bioinformatics 20, 696 (2019).
  * Karplus and McCammon (2002) M. Karplus and J. A. McCammon, Molecular dynamics simulations of biomolecules, Nat. Struct. Biol. 9, 646 (2002).
  * Münz  _et al._ (2012) M. Münz, J. Hein, and P. C. Biggin, The role of flexibility and conformational selection in the binding promiscuity of pdz domains, Plos Comput. Biol. 8, 1 (2012).
  * Pabon and Camacho (2017) N. A. Pabon and C. J. Camacho, Probing protein flexibility reveals a mechanism for selective promiscuity, eLife 6, e22889 (2017).
  * Kamenik  _et al._ (2021) A. S. Kamenik, I. Singh, P. Lak, T. E. Balius, K. R. Liedl, and B. K. Shoichet, Energy penalties enhance flexible receptor docking in a model cavity, P. Natl. Acad. Sci. Usa. 118, 10.1073/pnas.2106195118 (2021).
  * Petrović  _et al._ (2018) D. Petrović, V. A. Risso, S. C. L. Kamerlin, and J. M. Sanchez-Ruiz, Conformational dynamics and enzyme evolution, J. Roy. Soc. . Interface 15, 20180330 (2018).
  * Becker  _et al._ (2021) D. Becker, P. V. Bharatam, and H. Gohlke, F/g region rigidity is inversely correlated to substrate promiscuity of human cyp isoforms involved in metabolism, J. Chem. Inf. Model. 61, 4023 (2021).
  * Hill  _et al._ (2016) M. E. Hill, D. J. MacPherson, P. Wu, O. Julien, J. A. Wells, and J. A. Hardy, Reprogramming caspase-7 specificity by regio-specific mutations and selection provides alternate solutions for substrate recognition, Acs Chem. Biol. 11, 1603 (2016).
  * Eaton  _et al._ (1995) B. E. Eaton, L. Gold, and D. A. Zichi, Let’s get specific: The relationship between specificity and affinity, Chem. Biol. 2, 633 (1995).
  * Gao  _et al._ (2021) F. Gao, J. Glaser, and S. C. Glotzer, The role of complementary shape in protein dimerization, Soft Matter 17, 7376 (2021).
  * Greenspan (2010) N. S. Greenspan, Cohen’s conjecture, howard’s hypothesis, and ptashne’s ptruth: An exploration of the relationship between affinity and specificity, Trends Immunol. 31, 138 (2010).
  * Reyes and Kollman (2000) C. M. Reyes and P. A. Kollman, Structure and thermodynamics of rna-protein binding: Using molecular dynamics and free energy analyses to calculate the free energies of binding and conformational change, J. Mol. Biol. 297, 1145 (2000).
  * Kankainen  _et al._ (2004) M. Kankainen, T. Laitinen, and M. Peräkylä, Recognition of reactive high-energy conformations by shape complementarity and specific enzyme–substrate interactions in family 10 and 11 xylanases, Phys. Chem. Chem. Phys. 6, 5074 (2004).
  * Park  _et al._ (2016) J. Park, J. J. McDonald, R. C. Petter, and K. N. Houk, Molecular dynamics analysis of binding of kinase inhibitors to wt egfr and the t790m mutant, J. Chem. Theory Comput. 12, 2066 (2016).
  * Perona and Hadd (2012) J. J. Perona and A. Hadd, Structural diversity and protein engineering of the aminoacyl-trna synthetases, Biochemistry-us. 51, 8705 (2012).
  * Bar-Even  _et al._ (2015) A. Bar-Even, R. Milo, E. Noor, and D. S. Tawfik, The moderately efficient enzyme: Futile encounters and enzyme floppiness, Biochemistry-us. 54, 4969 (2015).
  * Eckmann and Tlusty (2021) J. P. Eckmann and T. Tlusty, Dimensional reduction in complex living systems: Where, why, and how, Bioessays 43, 2100062 (2021).
  * Ferreiro  _et al._ (2007) D. U. Ferreiro, J. A. Hegler, E. A. Komives, and P. G. Wolynes, Localizing frustration in native proteins and protein assemblies, P. Natl. Acad. Sci. Usa. 104, 19819 (2007).
  * Pande  _et al._ (1997) V. S. Pande, A. Y. Grosberg, and T. Tanaka, On the theory of folding kinetics for short proteins, Folding and Design 2, 109 (1997).
  * Ferreiro  _et al._ (2014) D. U. Ferreiro, E. A. Komives, and P. G. Wolynes, Frustration in biomolecules, Q. Rev. Biophys. 47, 285–363 (2014).
  * Consortium (2018) T. U. Consortium, Uniprot: A worldwide hub of protein knowledge, Nucleic Acids Res. 47, D506 (2018).
  * Bastolla and Demetrius (2005) U. Bastolla and L. Demetrius, Stability constraints and protein evolution: The role of chain length, composition and disulfide bonds, Protein Engineering, Design and Selection 18, 405 (2005).
  * Davidson  _et al._ (1995) A. R. Davidson, K. J. Lumb, and R. T. Sauer, Cooperatively folded proteins in random sequence libraries, Nat. Struct. Biol. 2, 856 (1995).
  * Maugh (1983) T. H. Maugh, Catalysts that break nature’s monopoly: Chiral complexes can approach the specificity of enzymes for synthesis of optically active compounds, and can act on a wider variety of substrates, Science 221, 351 (1983).
  * MacMillan (2008) D. W. C. MacMillan, The advent and development of organocatalysis, Nature 455, 304 (2008).
  * Henley  _et al._ (2020) M. J. Henley, B. M. Linhares, B. S. Morgan, T. Cierpicki, C. A. Fierke, and A. K. Mapp, Unexpected specificity within dynamic transcriptional protein–protein complexes, P. Natl. Acad. Sci. Usa. 117, 27346 (2020).
  * Savir and Tlusty (2008) Y. Savir and T. Tlusty, Optimal design of a molecular recognizer: Molecular recognition as a bayesian signal detection problem, Ieee J-stsp. 2, 390 (2008).
  * Echave (2008) J. Echave, Evolutionary divergence of protein structure: The linearly forced elastic network model, Chem. Phys. Lett. 457, 413 (2008).
  * Marcos and Echave (2020) M. L. Marcos and J. Echave, The variation among sites of protein structure divergence is shaped by mutation and scaled by selection, Current Research in Structural Biology 2, 156 (2020).
  * Tang and Kaneko (2021) Q. Y. Tang and K. Kaneko, Dynamics-evolution correspondence in protein structures, Phys. Rev. Lett. 127, 098103 (2021).
  * Byrd  _et al._ (1995) R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, A limited memory algorithm for bound constrained optimization, Siam J. Sci. Comput. 16, 1190 (1995).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Scikit-learn: Machine learning in Python, J. Mach. Learn. Res. 12, 2825 (2011).