# Binding-and-folding recognition of an intrinsically disordered protein using
online learning molecular dynamics

Pablo Herrera-Nieto  Adrià Pérez  Computational Science Laboratory,
Universitat Pompeu Fabra, Barcelona Biomedical Research Park (PRBB), C Dr.
Aiguader 88, 08003, Barcelona, Spain  Gianni De Fabritiis
<EMAIL_ADDRESS>Computational Science Laboratory, Universitat
Pompeu Fabra, Barcelona Biomedical Research Park (PRBB), C Dr. Aiguader 88,
08003, Barcelona, Spain

###### Abstract

Intrinsically disordered proteins participate in many biological processes by
folding upon binding with other proteins. However, coupled folding and binding
processes are not well understood from an atomistic point of view. One of the
main questions is whether folding occurs prior to or after binding. Here we
use a novel unbiased high-throughput adaptive sampling approach to reconstruct
the binding and folding between the disordered transactivation domain of c-Myb
and the KIX domain of the CREB-binding protein. The reconstructed long-term
dynamical process highlights the binding of a short stretch of amino acids on
c-Myb as a folded $\alpha$-helix. Leucine residues, specially Leu298 to
Leu302, establish initial native contacts that prime the binding and folding
of the rest of the peptide, with a mixture of conformational selection on the
N-terminal region with an induced fit of the C-terminal.

###### keywords:

Conformational selection, Coupled folding and binding, intrinsically
disordered proteins, molecular dynamics simulations, protein-protein
interactions

These authors contributed equally to this work Acellera Labs, C Dr Trueta 183,
08005, Barcelona, Spain These authors contributed equally to this work
Acellera Ltd, Devonshire House 582, HA7 1JS, United Kingdom Institució
Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluis Companys 23,
08010 Barcelona, Spain MD,MSM,IDP,UCB,NMR

##  1 Introduction

Intrinsically disordered proteins (IDPs) participate in many biological
functions despite lacking a stable tertiary structure 1. Initial clues for the
function of IDPs were revealed by structural studies 2, 3, showing that
proteins that were disordered in isolation became folded upon interacting with
their partners, opening to question how folding couples with binding.

Recently, molecular dynamics (MD) simulations have been successfully applied
to reconstruct biological dynamic events in problems such as protein-ligand 4
and protein-protein 5, 6 binding, as well as protein folding 7, 8. MD has also
been applied in the field of IDPs 9, 10, 11, 12. In particular, the Mdm2
protein and the disordered 12-residue N-terminal region of p53 were studied
using implicit solvent simulations9, parallel full-atom simulations totaling
11, biased free-energy-based sampling 10, and biased/unbiased simulations to
estimate kinetics on the second timescale 12. For another system, KIX-pKID, a
single event of binding 13 has been sampled at all-atom resolution.

The KIX—c-Myb binding-and-folding mechanism has been extensively studied
experimentally as an exemplar case of protein-IDP interaction 14, 15, 16, 17,
18, 19, 20. The KIX domain of the CREB-binding protein is a short 87-aa region
composed of three $\alpha$-helices (designated as $\alpha$-1, $\alpha$-2 and
$\alpha$-3, from N-terminal to C-terminal) forming a compact bundle 3. KIX
represents a paradigm of binding promiscuity: it binds to many IDPs, including
the proto-oncogene c-Myb 3 (Figure 1.a), with multiple binding conformations
14. However, the system composed by KIX—c-Myb remained outside of the scope of
all-atom molecular simulations due to the size of the IDP (it doubles the
length of p53) and the existence of multiple binding modes between them 14. In
particular, it is unclear whether the interaction takes place by
conformational selection, i.e. c-Myb needs to be folded before binding to its
partner, or by induced-fit, where binding not only happens independently of
c-Myb’s secondary structure but also triggers its folding, as shown for other
IDPs (KIX-pKID) 21, 13. Understanding these aspects has implications for the
druggability of disordered proteins. Another important factor is c-Myb’s high
helicity in isolation and the consequences it might exert on the final complex
structure, which features an extended $\alpha$-helical c-Myb bound to KIX.
Some reports support the induced-fit approach based on kinetics and
mutagenesis studies 16, 15, 20, while others advocate for a mixed mechanism
14; yet not a detailed model for the binding process is available.

In this paper, we take advantage of a novel algorithm that frames the MD
sampling problem from a reinforcement learning perspective (see 22 and
Methods) to reconstruct multiple binding modes between c-Myb and KIX. This
sampling algorithm was key for us to reconstruct the binding process, as
previous attempts over the years using other state-of-the-art adaptive
sampling methods 23, 24 were not successful, always failing to recover the NMR
bound structure. Results provide insights into the binding mechanism between
these two proteins, supporting a mixed model that combines both conformational
selection and induced fit.

Figure 1: Exploration performance. a) KIX—cMyb NMR structure. KIX domain is
shown as a white surface and ribbon and c-Myb bound to KIX as a red helix (PDB
code 1SB0). Exploration performance by b) Counts Adaptive ($\sim 480\ \mu s$)
, and AdaptiveBandit ($\sim 450\ \mu s$) is shown by plotting the mean RMSD
(on the x-axis) and standard deviation (on the y-axis) for each of the MSM’s
microstates, color-mapped accordingly to their macrostate assignment. The
dashed square indicates the bound zone, placed in the region corresponding to
low mean RMSD and standard deviation.

##  2 Results & Discussion

### Adaptive sampling the KIX—c-Myb binding-and-folding process

Simulations to reconstruct the KIX—c-Myb binding mode were performed following
an adaptive sampling strategy. In adaptive sampling, successive rounds of
simulations are performed in an iterative step-wise manner, where an
acquisition function over the currently sampled conformations is defined. We
compare two of the acquisition functions used for KIX—c-Myb simulations: a
count-based one and another one inspired by reinforcement learning, part of
the novel AdaptiveBandit method 22.

The new AdaptiveBandit method is framed into a simplified reinforcement
learning problem, the multi-armed bandit problem (see Methods). We use the
upper confidence bound (UCB) algorithm 25 to optimize an action-picking policy
in order to maximize future rewards, optimally balancing the exploration of
new higher rewarding actions with the exploitation of the most known rewarding
ones. The reward function, which associates the action with the reward given
by the system, defines what we want to optimize. In this work, we choose the
reward to be minus the free energy of each configuration visited in the
trajectory spawn from a given action (see Eq.2 in Methods), where the free
energy of a conformation is given by the corresponding Markov state model
(MSM) microstate computed with the data available at the current sampling
epoch.

Standard low counts adaptive sampling 23 (hereby named Counts Adaptive) can be
shown to be optimal in pure exploration conditions 24. Counts are computed
over clusters of conformations; this method is, however, noisy as clusters can
be poorly populated. Therefore, in the implementation available in 24, counts
are computed over a smaller subset by grouping clusters (microstates) into
macrostates, constructing a Markov State Model (MSM) 26 with the available
data at each round. The acquisition function is given by proportionally
choosing macrostates as $1/c$, where $c$ represents macrostate counts, and by
randomly selecting conformations within them.

A comparison between Counts Adaptive and AdaptiveBandit is provided in Figure
1.b. The batch based on Counts Adaptive (48 epochs) failed to connect
microstates similar to the NMR structure in over $\sim 480\ \mu s$, reaching
at best an RMSD around 7 Å. For us, it was impossible to build an MSM with the
bound state with previous methods, and novel approaches were needed to
reconstruct the binding-and-folding process between KIX and c-Myb
successfully. AdaptiveBandit provides converged estimates of kinetics and
thermodynamics after just $150\ \mu s$ of sampling (Supplementary Figure 5).

Figure 2: KIX & c-Myb binding model. a) States distribution across the TICA
space: microstates are represented as dots and are colored following their
macrostate assignment. b) Representative structures: PDB structure 1SB0 is
depicted with KIX as a gray surface, c-Myb bound to the primary interface as a
yellow ribbon, and c-Myb bound to the secondary interface as a blue ribbon.
c-Myb backbones for 30 representatives MD structures of bound and secondary
states are displayed with blurry yellow and blue clouds, respectively. c)
Macrostate contact fingerprint: profile of contacts established between c-Myb
and KIX in each macrostate in at least 50% of the structures. Blue color
represents contacts present in the state but not in the original NMR
structure; green indicates original NMR contacts not found in the MSM state;
and yellow squares represent contact matches, found in both NMR and MD
structures. d) Macrostate cMyb helicity: helicity fraction per residue of
c-Myb in each macrostate. Helicity for the cMyb peptide alone is depicted in
grey in each plot for comparison.

### Identification of the bound state

The full data set of the AdaptiveBandit run accounted for a total simulation
time of $\sim$450 $\mu$s, split across 40 epochs, and was the one used to
study the molecular features of KIX—c-Myb binding-and-folding. For the
analysis, a slightly different MSM was built based on all-pair
$C_{\alpha}+C_{\beta}$ distances between KIX and c-Myb, self distances between
$C_{\alpha}$ of c-Myb, secondary structure of c-Myb and RMSD to the NMR bound
conformation (PDB ID: 1SB0). The MSM defines three kinetically similar sets of
conformations, referred as macrostates (Figure 2.a and Supplementary Figure
1.b): a highly populated state with an heterogeneous mixture of conformations
(unbound), a well defined c-Myb bound state (bound) and, finally, a secondary
bound state (secondary). Representative structures of all states can be found
in Figure 2.b. The bound macrostate contains structures with a minimum RMSD of
approximately 3 Å with respect to the NMR structure. Complete binding and
folding trajectory videos can be found in Supplementary Table 1, with
reconstructed trajectories from different epochs containing unique paths to
the bound state (with RMSD $<4$ Å with respect to the bound NMR structure)

The bound state identifies the primary cMyb bound pose in the hydrophobic
groove between $\alpha$-1 and $\alpha$-3) of KIX. On average, it shares 36% of
the fraction of native inter-molecular contacts ($Q_{int}$) with the original
NMR structure, as shown in Figure 2.c. These contacts mainly involve the
interaction of c-Myb residues Leu298 and Leu302 with residues across the
primary binding interface: Leu302 contacts Leu603, Leu653, and specially
Leu607 of KIX, which is buried down in the pocket, whereas Leu298 establishes
additional native contacts with Ala610, Ile657, and Tyr658. $Q_{int}$ reaches
up to 80% in those microstates exhibiting the tightest bound conformations,
and, in addition to the leucine binding, they feature most of the contacts
between the C-terminal half of c-Myb and KIX, which are not that prevalent
across the bound macrostate (Supplementary Figure 2). The main contacts
missing account for the electrostatic interactions established between Arg294
and the region on $\alpha$-3. There are some conformations where these
interactions occur, but their prevalence in those microstates is less than
50%.

The secondary structure profile for MD-derived states matches the experimental
description of c-Myb 14, 19, as shown in Supplementary Figure 3: the 25
residues are separated in two halves by residues Met303 and Ser304. The
N-terminal half shows a high helical tendency, around 20-30% for residues in
positions 297 to 302 with c-Myb in isolation, being maximal in bound states.
Experimentally, this N-terminal half in isolation reaches even higher helicity
levels ($\sim$70%) when using an extended construct of c-Myb 14. On the other
hand, the C-terminal section exhibits low helical propensity when in
isolation, and increases when bound to KIX. The full helix conformation only
appears in those microstates with the tightest bound conformations.

Figure 3: Complete c-Myb binding process to KIX domain. Main pathways leading
from Bulk (macrostate 14) to the Bound state (macrostate 12). Fluxes are shown
as percentages near the arrows. Only those fluxes higher than 5% are shown.
The arrow thickness is proportional to the flux percentage. Straight arrows
indicate the maximum flux path, while dashed arrows show other fluxes. Each
macrostate structure shows KIX as the white surface and ribbons and c-Myb as
the orange ribbon and tubes. For each macrostate, 25 conformations are shown
as thin tubes with one structure highlighted as a ribbon structure that
includes the side chain of Leu302, colored by atom element. Additionally, for
macrostates 7 and 12, the reference NMR c-Myb structure is shown as a
transparent red ribbon for comparison.

### Secondary binding mode

The existence of alternative binding poses between c-Myb and KIX has also been
reported 14. The MSM shows the presence of a secondary binding mode (referred
to as secondary), occupying a novel interface, located between $\alpha$-1 and
$\alpha$-2 (Figure 1.b and Supplementary Figure 9). The interaction of the
secondary state resembles the bound binding mode: the N-terminal half is
folded in the typical $\alpha$-helix, while the C-terminal section remains
mostly unstructured. The presence of a native contact in this secondary
binding mode is due to the penetration of Leu302, located close to Leu603’s
backbone in KIX, rather than by side-chain proximity. Leu298 and Leu302 of
c-Myb are deeply buried in a hydrophobic pocket composed of residues Val604,
Val608, Leu620 (found in the G2 helix, which connects $\alpha$-1 and
$\alpha$-2) and Val629. Kinetically, there is a 10-fold difference in the mean
first passage time for binding between both sites — $(9.96\pm 3.57)\cdot
10^{3}$ ns for binding to bound site and $(1.05\pm 0.46)\cdot 10^{5}$ ns for
the secondary site — that may account for the preferential binding of c-Myb to
the primary interface.

### Model validation

To validate the model, we compared the kinetic parameters derived from it with
available information 17. Experimental values from Shammas et al. were
calculated at temperatures ranging from 278 to 298 K, while simulations were
executed at physiological temperature (310 K). $k_{on}$ values display a
temperature independent tendency, whereas for temperature dependent variables
$k_{off}$ and $k_{d}$ values had to be extrapolated to 310K (Supplementary
Figure 4). Hence, reference values for $k_{off}$ and free energy (obtained
from $k_{d}$) resulted in 866 $s^{-1}$ and -6.81 $kcal\ mol^{-1}$
respectively.

Due to the size of the peptide compared to the solvation box, it is hard for
the MSM to automatically define the correct bulk state. We, therefore,
manually defined a bulk state that contains conformations where the distance
between KIX and cMyb is maximized. The bulk state was defined by taking those
microstates where the minimum distance between KIX and cMyb is higher than a
threshold. Consequently, some kinetic and thermodynamic estimates have a
dependency on such distance threshold (Supplementary Figure 6a) as this
affects the definition of the bulk state. However, the computed $k_{off}$ and
free energy estimates are practically stable after a minimum separation
distance of just 4 Å.

The obtained MSM estimations of $k_{on}$ go between $(2.72)\cdot
10^{7}M^{-1}s^{-1}$ and $(3.65)\cdot 10^{7}M^{-1}s^{-1}$, in agreement with
the experimental value $(2.2\pm 0.2)\cdot 10^{7}M^{-1}s^{-1}$ 17. $k_{off}$
estimates range from $3.50\cdot 10^{3}\ s^{-1}$ to $21.70\cdot 10^{3}\
s^{-1}$, overestimating the extrapolated experimental value by an order of
magnitude. Free energy estimates range between $-7.35\ kcal\ mol^{-1}$ and
$-6.27\ kcal\ mol^{-1}$ depending on the choices of the analysis parameters,
with the extrapolated experimental value inside this interval.

We further verified the reproducibility of kinetic and thermodynamic
measurements to ensure model convergence by building multiple MSMs using
incrementally more trajectories. Convergence is reached at $150\ \mu s$ on all
the previous estimates (Supplementary Figure 5). The discrepancy between the
experimental reference and computed $k_{off}$ values translates into a faster
dissociation in our model. We also verified if this was due to normal
discretization errors in the MSM projection or to the fact that our
simulations did not obtain a complete bound conformation between KIX and cMyb.
In order to test this hypothesis, additional long trajectories (8 replicas of
2 $\mu$s each) were run starting from bound NMR and MD-derived conformations.
We constructed an MSM using both simulation datasets. However, the free energy
estimations are only marginally improved (Supplementary Figure 6b). Thus, we
concluded that the additional bound simulations do not add additional
information and we restrict the analysis to just the AdaptiveBandit set of
simulations as this is the most general case where no NMR information is
available.

### Binding follows both induced-fit and conformational selection

In order to gain additional structural insight of the binding process, we
constructed an MSM with a higher number of macrostates, using the same lag
time. We used transition path theory 27, 28 to calculate fluxes leading from
the purely bulk state to the bound conformations. Out of the 15 macrostates of
the new MSM, only a reduced set of 6 is sufficient to explain binding to the
primary interface. The other macrostates describe either the secondary binding
mode or other unstable interactions between KIX and cMyb. The network
generated by the flux interchanges between macrostates (Figure 3) separates
the binding and folding process into three events: the establishment of the
initial contacts, binding and folding of the N-term section of cMyb and
finally, binding and folding of the C-term section of cMyb. The first step of
the binding process features the first native contacts found across the
KIX—c-Myb binding pathway, which involves residues Leu302 of c-Myb. The role
of Leu302 as the main driving force for the interaction has already been
described 3 and is due in part to the kink in the helix created by neighbors
residues Met303 and Ser304, which exposes Leu302 allowing for a deep
penetration inside the binding pocket. Besides, of the KIX residues contacted
at this stage is Leu603, which is one of the most exposed residues in the
hydrophobic pocket later occupied by Leu302.

To determine if cMyb folding precedes or follows binding at this step, we
looked at the flux passing through different cMyb conformations in macrostate
13 (Supplementary Figure 8). We see that almost half of the flux goes through
N-term helical conformations, suggesting that the presence of helix in
residues 297 to 302 facilitates these first binding step, following a
conformational selection mechanism. There is also a considerable flux going
through unfolded conformations, meaning Leu302 also binds through induced fit.

The second step of the binding process goes from macrostate 13 to macrostate
11, where several contacts are formed across the N-term of cMyb. The last
step, which goes from macrostate 11 to 12 (the bound state), involves forming
the last contacts on the C-term and completely folding cMyb. Here, we also
looked at the flux passing through macrostate 11 (Supplementary Figure 8c) to
discern between conformational selection and induced fit on the C-term folding
and binding mechanism. Here, all the flux goes through conformations where the
C-term is unfolded, meaning that only when the C-term native contacts start to
happen we see a complete cMyb folding to an alpha helix, following a clear
induced-fit mechanism.

The overall mechanism works as an induced-fit binding and folding, but we can
see a mixed mechanism during the first binding steps, where both
conformational selection and induced-fit seem to take a part in facilitating
the first contacts through cMyb’s Leu302. In summary, initial steps can be
greatly benefited from pre-folded helical structures of c-Myb (Figure 3),
although binding before folding is also observed. The binding of helical
conformations dominates the initial steps of the interaction, but for the
interaction of the C-terminal tail, folding follows binding. No limiting steps
in the binding process are observed; hence no possible transition states can
be defined, as pointed out by experimental reports 17.

##  3 Conclusions

The analysis presented here provides a detailed molecular description of
binding of c-Myb to the primary interface of KIX, summarized as a two-step
process, where initially the N-terminal region of c-Myb binds with a preferred
helical conformation, allowing the formation of native contacts and, in the
last step, folding and binding of the C-terminal. Study of the fluxes derived
from the MSM shows the relevance of residue Leu302, not only in the final
bound structure but also as the responsible for establishing the first
contacts and serving as an anchoring point between c-Myb and KIX.

The model describes an overall induced-fit binding mechanism, as the complete
folding of cMyb is only observed when native contacts have been formed.
Conformational selection would only affect the first binding stage on residues
298 to 302 and not the whole length of the peptide, whereas the latter stages
of binding follow an induced-fit mechanism.

Overall, our results provide a detailed mechanistic model for the binding of
c-Myb to the primary interface of KIX, as well as showing the interaction with
a secondary binding site, by using unbiased full-atom MD simulations and MSM
analysis. The novel MD sampling approach used in this work, AdaptiveBandit,
had a crucial role in resolving this type of folding and binding process. The
method is implemented and available in the HTMD python package 24. However,
more algorithms can be derived within the same bandit framework. While here we
choose the reward to be minus the free energy, other choices could optimize
different costs, for example, improving the precision of the off-rate or
optimizing sampling in the context of structure prediction.

##  4 Methods

### Molecular dynamics simulations

In order to generate initial conformations for c-Myb (residues 291 to 315), we
ran multiple parallel simulations. The peptide was solvated in a cubic water
box of $64$ Å side with a NaCl concentration of 0.05 $M$. First, the peptide
was simulated at 500 $K$ for 120 ns to unfold the initial structure. Then, 200
systems were built by placing one random unstructured c-Myb conformation in
conjunction with KIX in opposite corners of a 64 Å side cubic water box with a
NaCl concentration of 0.05 $M$, resulting in a final protein concentration of
$\sim$3.2 mM.

All systems were built using HTMD 24 and simulated with ACEMD 29, CHARMM22*
force field 30 and TIP3P water model 31. A Langevin integrator was used with a
damping constant of 0.1 ps-1. The integration time step was set to 4 fs, with
heavy hydrogen atoms (scaled up to four times the hydrogen mass) and holonomic
constraints on all hydrogen-heavy atom bond terms. Electrostatics were
computed using PME with a cutoff distance of 9 Å and grid spacing of 1 Å.
After energy minimization, equilibration for all systems was done in an NPT
ensemble at 303 K, 1 atm, with heavy atoms constrained at 1 $kcal\ mol^{-1}$
Å2. Energy minimization was run for 500 steps and equilibrated for 2 ns.

Production runs of 250 $ns$ were performed at 310 $K$ using the distributed
computing project GPUGrid 32, following an adaptive sampling strategy. The
final data set included 1,809 trajectories of 250 ns, resulting in an
aggregated simulation time of $\sim$450 $\mu$s. Additionally, a set of long MD
runs were performed starting from bound structures. Four models of the NMR-
determined structure and four random bound conformations were selected and
equilibrated as previously described. A total of 8 long trajectories of $\sim
2\ \mu s$ each were generated.

### Markov state model analysis

The projected space used for building the MSM included four different
featurizations: all pair Cα \+ Cβ atoms distances between KIX and c-Myb to
account for the interaction between the two proteins, self-distances between
every Cα of c-Myb and its secondary structure, to monitor its conformation,
and finally, RMSD to the bound structure. TICA was used at a lag time
$\tau=20$ ns (implied timescales are shown in Supplementary Figure 1.a) for
both the distance features and the secondary structure features, taking the 4
most relevant components from the distance features (both inter-distances and
cMyb self-distances) and the 3 most relevant components from the secondary
structure features.

The 8-dimensional projected data was discretized into 2,000 clusters using the
mini-batch k-means algorithm 33. The microstates defined in the MSM were
coarse-grained into larger meta-stable macrostates by using PCCA++ 34. For the
estimation of kinetic values, the original MSM was modified by creating an
additional macrostate, considered as the bulk state for all subsequent
calculations to obtain the kinetics of binding. The bulk state was created by
taking those microstates where the minimum distance between KIX and cMyb was
higher than a threshold. Error in kinetic measures was estimated by creating
50 independent MSMs using a random set containing 80% of the simulation data.

To obtain the kinetic pathway of binding and folding, we increased the number
of macrostates in the MSM using PCCA++ again. Fluxes between macros were
estimated using transition path theory 27, 28. For the intra-macrostate flux
analysis, we computed the mean helicity of cMyb for each microstate in it and
clustered them into 4 main states which describe the peptide’s grade of helix
formation. All analysis were performed with HTMD 24.

### AdaptiveBandit sampling

The multi-armed bandit problem is defined by
$\langle\mathcal{A},\mathcal{R},\gamma\rangle$, where an action
$a_{t}\in\mathcal{A}$ and $\mathcal{R}^{a}$ is a (stochastic) reward function.
We choose $\gamma=0$ for totally discounted rewards. The optimal policy
$\pi_{a}\sim\mathbb{P}[a]$ selects actions $a_{t}$ in order to maximize the
cumulative future rewards. The construction of an optimal selection strategy
requires handling the exploration-exploitation problem. AdaptiveBandit relies
on the UCB1 algorithm 25, defining an upper confidence bound for each action-
value estimate based on the number of times an action has been picked and the
total amount of actions taken

$a_{t}=\operatorname*{argmax}_{a\in\mathcal{A}}\left[{Q_{t}(a)+c\sqrt{\frac{\ln{t}}{N_{t}(a)}}}\right],$
(1)

where $t$ denotes the total number of actions taken,
$Q_{t}(a)=\mathbb{E}_{\pi}[r|a]$ is the action-value estimation, $N_{t}(a)$ is
the number of times action $a$ has been selected (prior to time $t$) and $c$
is a constant controlling the degree of exploration. As for the reward
definition, there are different choices depending on the objective, e.g. here,
the interest is sampling the bound metastable state, hence, we rewarded
actions based on the stability of conformations using MSM estimations of the
free energy for each state

$\mathcal{R}_{a}=<k_{B}T\log(\mu(x))>_{(a,x_{1},\dots,x_{\tau})},$ (2)

where $\mu(x)$ is the equilibrium distribution estimated by the MSM with the
currently available data and the average is performed over the frames in the
trajectory starting from $a$. AdaptiveBandit uses the MSM discretized
conformational space to define the action set and at each round acquires a
random conformation from the selected states to respawn new simulations. A
more formal description of the bandit framework and AdaptiveBandit in the
context of adaptive sampling as well as analysis in simpler, analytical
potentials are available at 22. The AdaptiveBandit sampling algorithm is made
available in the HTMD 24 Python package.

### Adaptive Sampling parameters

For both the AdaptiveBandit and the count Adaptive runs, the construction of
MSMs at each epoch was done using the residue-residue contacts between KIX and
c-Myb measured as the minimum contacts between residues at a threshold of 5 Å,
and the backbone dihedral angles of c-Myb. Time independent component analysis
(TICA) 35 was used for dimensionality reduction using a lag time of $\tau=20$
frames and keeping the 3 first dimensions, which were later clustered with a
k-centers algorithm. AdaptiveBandit was performed during 40 epochs with a $c$
value of 0.01.

The authors thank Kresten Lindorff-Larsen and Frank Noé for their critical
reading of the manuscript. The authors also thank volunteers at GPUGRID.net
for contributing computational resources and Acellera for funding. G.D.F.
acknowledges support from MINECO (Unidad de Excelencia María de Maeztu
CEX2018-000782-M and BIO2017-82628-P) and FEDER. This project received funding
from the European Union’s Horizon 2020 Research and Innovation Programme under
Grant Agreement No. 823712 (CompBioMed2 Project).

##  5 Author Contributions

P.H.N. and A.P. generated and analyzed the data; P.H.N., A.P. and G.D.F. wrote
the paper; G.D.F. designed research.

Supporting information is available free of charge.

  * •

Supporting information: Contains all supporting figures as well as the
supporting table containing trajectory videos capturing the folding and
binding of c-Myb with KIX domain.

## References

  * Dyson and Wright 2005 Dyson, H. J.; Wright, P. E. Intrinsically unstructured proteins and their functions. _Nature reviews Molecular cell biology_ 2005, _6_ , 197
  * Kussie et al. 1996 Kussie, P. H.; Gorina, S.; Marechal, V.; Elenbaas, B.; Moreau, J.; Levine, A. J.; Pavletich, N. P. Structure of the MDM2 oncoprotein bound to the p53 tumor suppressor transactivation domain. _Science_ 1996, _274_ , 948–953
  * Zor et al. 2004 Zor, T.; De Guzman, R. N.; Dyson, H. J.; Wright, P. E. Solution structure of the KIX domain of CBP bound to the transactivation domain of c-Myb. _Journal of molecular biology_ 2004, _337_ , 521–534
  * Buch et al. 2011 Buch, I.; Giorgino, T.; De Fabritiis, G. Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations. _Proceedings of the National Academy of Sciences_ 2011, _108_ , 10184–10189
  * Plattner et al. 2017 Plattner, N.; Doerr, S.; De Fabritiis, G.; Noé, F. Complete protein–protein association kinetics in atomic detail revealed by molecular dynamics simulations and Markov modelling. _Nature chemistry_ 2017, _9_ , 1005
  * Borgia et al. 2018 Borgia, A.; Borgia, M. B.; Bugge, K.; Kissling, V. M.; Heidarsson, P. O.; Fernandes, C. B.; Sottini, A.; Soranno, A.; Buholzer, K. J.; Nettels, D., et al. Extreme disorder in an ultrahigh-affinity protein complex. _Nature_ 2018, _555_ , 61
  * Lindorff-Larsen et al. 2011 Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. How fast-folding proteins fold. _Science_ 2011, _334_ , 517–520
  * Piana et al. 2013 Piana, S.; Lindorff-Larsen, K.; Shaw, D. E. Atomistic description of the folding of a dimeric protein. _The Journal of Physical Chemistry B_ 2013, _117_ , 12935–12942
  * Zwier et al. 2016 Zwier, M. C.; Pratt, A. J.; Adelman, J. L.; Kaus, J. W.; Zuckerman, D. M.; Chong, L. T. Efficient atomistic simulation of pathways and calculation of rate constants for a protein–peptide binding process: application to the MDM2 protein and an intrinsically disordered p53 peptide. _The journal of physical chemistry letters_ 2016, _7_ , 3440–3445
  * Morrone et al. 2017 Morrone, J. A.; Perez, A.; MacCallum, J.; Dill, K. A. Computed binding of peptides to proteins with MELD-accelerated molecular dynamics. _Journal of chemical theory and computation_ 2017, _13_ , 870–876
  * Zhou et al. 2017 Zhou, G.; Pantelopulos, G. A.; Mukherjee, S.; Voelz, V. A. Bridging microscopic and macroscopic mechanisms of p53-MDM2 binding with kinetic network models. _Biophysical journal_ 2017, _113_ , 785–793
  * Paul et al. 2017 Paul, F.; Wehmeyer, C.; Abualrous, E. T.; Wu, H.; Crabtree, M. D.; Schöneberg, J.; Clarke, J.; Freund, C.; Weikl, T. R.; Noé, F. Protein-peptide association kinetics beyond the seconds timescale from atomistic simulations. _Nature communications_ 2017, _8_ , 1095
  * Chong et al. 2019 Chong, S.-H.; Im, H.; Ham, S. Explicit Characterization of the Free Energy Landscape of pKID–KIX Coupled Folding and Binding. _ACS Central Science_ 2019, _5_ , 1342–1351
  * Arai et al. 2015 Arai, M.; Sugase, K.; Dyson, H. J.; Wright, P. E. Conformational propensities of intrinsically disordered proteins influence the mechanism of binding and folding. _Proceedings of the National Academy of Sciences_ 2015, _112_ , 9614–9619
  * Giri et al. 2013 Giri, R.; Morrone, A.; Toto, A.; Brunori, M.; Gianni, S. Structure of the transition state for the binding of c-Myb and KIX highlights an unexpected order for a disordered system. _Proceedings of the National Academy of Sciences_ 2013, _110_ , 14942–14947
  * Gianni et al. 2012 Gianni, S.; Morrone, A.; Giri, R.; Brunori, M. A folding-after-binding mechanism describes the recognition between the transactivation domain of c-Myb and the KIX domain of the CREB-binding protein. _Biochemical and biophysical research communications_ 2012, _428_ , 205–209
  * Shammas et al. 2013 Shammas, S. L.; Travis, A. J.; Clarke, J. Remarkably fast coupled folding and binding of the intrinsically disordered transactivation domain of cMyb to CBP KIX. _The journal of physical chemistry B_ 2013, _117_ , 13346–13356
  * Toto et al. 2016 Toto, A.; Camilloni, C.; Giri, R.; Brunori, M.; Vendruscolo, M.; Gianni, S. Molecular recognition by templated folding of an intrinsically disordered protein. _Scientific reports_ 2016, _6_ , 21994
  * Poosapati et al. 2018 Poosapati, A.; Gregory, E.; Borcherds, W. M.; Chemes, L. B.; Daughdrill, G. W. Uncoupling the folding and binding of an intrinsically disordered protein. _Journal of molecular biology_ 2018, _430_ , 2389–2402
  * Shammas et al. 2014 Shammas, S. L.; Travis, A. J.; Clarke, J. Allostery within a transcription coactivator is predominantly mediated through dissociation rate constants. _Proceedings of the National Academy of Sciences_ 2014, _111_ , 12055–12060
  * Sugase et al. 2007 Sugase, K.; Dyson, H. J.; Wright, P. E. Mechanism of coupled folding and binding of an intrinsically disordered protein. _Nature_ 2007, _447_ , 1021
  * Pérez et al. 2020 Pérez, A.; Herrera-Nieto, P.; Doerr, S.; De Fabritiis, G. AdaptiveBandit: a multi-armed bandit framework for adaptive sampling in molecular simulations. _Journal of Chemical Theory and Computation_ 2020, _16_ , 4685–4693
  * Doerr and De Fabritiis 2014 Doerr, S.; De Fabritiis, G. On-the-fly learning and sampling of ligand binding by high-throughput molecular simulations. _Journal of chemical theory and computation_ 2014, _10_ , 2064–2069
  * Doerr et al. 2016 Doerr, S.; Harvey, M.; Noé, F.; De Fabritiis, G. HTMD: high-throughput molecular dynamics for molecular discovery. _Journal of chemical theory and computation_ 2016, _12_ , 1845–1852
  * Auer 2002 Auer, P. Using confidence bounds for exploitation-exploration trade-offs. _Journal of Machine Learning Research_ 2002, _3_ , 397–422
  * Prinz et al. 2011 Prinz, J.-H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J. D.; Schütte, C.; Noé, F. Markov models of molecular kinetics: Generation and validation. _The Journal of chemical physics_ 2011, _134_ , 174105
  * Weinan and Vanden-Eijnden 2006 Weinan, E.; Vanden-Eijnden, E. Towards a theory of transition paths. _Journal of statistical physics_ 2006, _123_ , 503
  * Noé et al. 2009 Noé, F.; Schütte, C.; Vanden-Eijnden, E.; Reich, L.; Weikl, T. R. Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations. _Proceedings of the National Academy of Sciences_ 2009, _106_ , 19011–19016
  * Harvey et al. 2009 Harvey, M. J.; Giupponi, G.; Fabritiis, G. D. ACEMD: accelerating biomolecular dynamics in the microsecond time scale. _Journal of chemical theory and computation_ 2009, _5_ , 1632–1639
  * Piana et al. 2011 Piana, S.; Lindorff-Larsen, K.; Shaw, D. E. How robust are protein folding simulations with respect to force field parameterization? _Biophysical journal_ 2011, _100_ , L47–L49
  * Jorgensen et al. 1983 Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. _The Journal of chemical physics_ 1983, _79_ , 926–935
  * Buch et al. 2010 Buch, I.; Harvey, M. J.; Giorgino, T.; Anderson, D. P.; De Fabritiis, G. High-throughput all-atom molecular dynamics simulations using distributed computing. _Journal of chemical information and modeling_ 2010, _50_ , 397–403
  * Pedregosa et al. 2011 Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V., et al. Scikit-learn: Machine learning in Python. _Journal of machine learning research_ 2011, _12_ , 2825–2830
  * Röblitz and Weber 2013 Röblitz, S.; Weber, M. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. _Advances in Data Analysis and Classification_ 2013, _7_ , 147–179
  * Pérez-Hernández et al. 2013 Pérez-Hernández, G.; Paul, F.; Giorgino, T.; De Fabritiis, G.; Noé, F. Identification of slow molecular order parameters for Markov model construction. _The Journal of chemical physics_ 2013, _139_ , 07B604_1