# Exploring PROTAC cooperativity with coarse-grained alchemical methods

Huanghao Mai  [<EMAIL_ADDRESS>Matthew H. Zimmer  [  Thomas F. Miller III
[

###### Abstract

Proteolysis targeting chimera (PROTAC) is a novel drug modality that
facilitates the degradation of a target protein by inducing proximity with an
E3 ligase. In this work, we present a new computational framework to model
the cooperativity between PROTAC-E3 binding and PROTAC-target binding
principally through protein-protein interactions (PPIs) induced by the PROTAC.
Due to the scarcity and low resolution of experimental measurements, the
physical and chemical drivers of these non-native PPIs remain to be
elucidated. We develop a coarse-grained (CG) approach to model interactions in
the target-PROTAC-E3 complexes, which enables converged thermodynamic
estimations using alchemical free energy calculation methods despite an
unconventional scale of perturbations. With minimal parameterization, we
successfully capture fundamental principles of cooperativity, including the
optimality of intermediate PROTAC linker lengths that originates from
configurational entropy. We qualitatively characterize the dependency of
cooperativity on PROTAC linker lengths and protein charges and shapes. Minimal
inclusion of sequence- and conformation-specific features in our current
forcefield, however, limits quantitative modeling to reproduce experimental
measurements, but further development of the CG model may allow for efficient
computational screening to optimize PROTAC cooperativity.

Caltech] Division of Chemistry and Chemical Engineering, California Institute
of Technology Caltech] Division of Chemistry and Chemical Engineering,
California Institute of Technology Caltech] Division of Chemistry and Chemical
Engineering, California Institute of Technology

##  1 Introduction

Proteolysis targeting chimera (PROTAC) has emerged as a promising drug
modality that elicits protein degradation by hijacking the ubiquitin-
proteasome system (UPS), a major regulatory component of cells. In the UPS
pathway, E3 ligases transfer ubiquitins onto aberrant proteins to mark them
for degradation by proteasomes. A PROTAC molecule exploits this pathway with
two binding moieties that tether the target protein and an E3 ligase together.
The tethered target protein thus becomes a neo-substrate of the E3 ligase and
is subsequently ubiquitinated for proteasomal degradation. PROTACs require a
lower dose than conventional small-molecule inhibitors because of their
catalytic nature and they have the potential to target the undruggable
proteome 1, 2. Since the first proof-of-concept in 2001 3, the number of
proteins successfully degraded by PROTACs has grown rapidly, and examples of
such proteins include kinases and gene regulators that are implicated in
cancer. As of 2021, at least 13 PROTACs are in or approaching clinical trials
4.

Despite increasing applications, there is a lack of guidance on designing
PROTACs due to the unique mode of action 5, 6, 7. In particular, a critical
step in the degradation process is the formation of the ternary complex of
target-PROTAC-E3. The ternary complex involves molecular interactions beyond
the binary bindings between the two warheads of a PROTAC and the two proteins.
The selectivity 8, 9, 10 and stability 11, 12, 13, 14 of the ternary complex
can both be improved through favorable  protein-protein interactions (PPIs)
between the target protein and the E3 ligase. For certain targets, the
degradation outcome can be very different depending on whether cereblon (CRBN)
or von Hippel-Lindau (VHL), the two most heavily used E3 ligases, more
efficiently and selectively form a productive complex with the target 15, 16,
11, 17. As more warheads for E3 ligases are designed 18, 19, 20, 21, choosing
which of the more than $600$ E3 ligases in humans 22 optimally interact with
the target protein will become important 23, 24. While PPIs depend on the
sequences and the structures of the proteins, PROTACs can also modulate the
PPIs by restricting the distance and relative orientation between the target
and the E3 ligase, effectively changing the entropic component of PPIs.

Because of this three-body interplay and the transient nature of the ternary
complex, a complete characterization of the PPIs as a function of the PROTAC,
the target protein, and the E3 ligase is intractable. A few proteomics studies
16, 25, 17 on kinase degradation have used PROTACs with promiscuous warheads
such that the PROTAC-induced PPIs differentially affect the degradation
outcome of hundreds of proteins. These studies reported the fold change of
protein abundance due to PROTAC treatment, but analysis can be complicated by
secondary interactions 24 and numerous other factors such as the permeability
of the PROTAC, half-lives of the target proteins, cellular localization, and
reactions downstream of ternary complex formation 26. Other studies 8, 27, 9,
28, 29, 30 have focused on specific target-E3 pairs and examined the effect of
changing PROTAC properties such as the linker length. They measured the
difference in the strength of PROTACs binding to the target or the E3 ligase
due to the presence of the other protein. This difference, termed binding
cooperativity, reflects the strength of PROTAC-mediated PPIs. However, few
generalizable patterns have emerged and systematic experimental
characterizations remain scarce.

Computational modeling based on docking or atomistic molecular dynamics (MD)
has complemented experimental work 9, 29 and displayed promising future
prospects, but there are several limitations to current methodologies.
Although standard docking protocols don’t handle three-body problems, several
workflows have been adapted ad hoc for PROTAC 31, 32, 33, 34, 35. Docking
studies rank ternary complex conformations by scoring functions biased for
naturally evolved PPIs and benchmark against the few crystal structures of
PROTAC-induced ternary complexes 36, 37, 38. The results can be inaccurate as
PROTAC-induced PPIs are non-native and exhibit plasticity 9, 39. In contrast,
atomistic MD is physically grounded to capture non-native PPIs. However, the
size of the ternary complex modeled at an atomistic resolution significantly
limits the timescale of simulations, such that naively simulating PPIs can be
prohibitively slow. Sophisticated enhanced sampling techniques and distributed
computing are needed to sample an ensemble of low-energy conformations that
are consistent with experimental data 40. Due to the difficulties in modeling
the ternary complex, direct calculation of the binding cooperativities was not
attempted until two recent studies 41, 42 that explored the molecular
mechanics with the generalized Born and surface area continuum solvation
(MM/GBSA).

Here, we seek an orthogonal approach that combines coarse-grained MD (CGMD)
and alchemical free energy calculation methods to study PROTAC
cooperativities. On the spectrum of computational tools, docking and atomistic
MD are positioned at the empirical and first-principle ends respectively, and
finding a compromise in the middle of this spectrum is a promising direction.
Compared to atomistic modeling, coarse-graining reduces the effective size of
the model and smoothens the energy surface, enabling simulations at a much
longer timescale necessary for the PROTAC-mediated complexes. While CGMD may
struggle to recapitulate the molecular basis of lock-and-key bindings, such
strong and specific interactions are less imperative in non-native PPIs
induced by PROTACs. Moreover, PROTAC binding reduces the ways proteins can
interact with each other, differentiating and simplifying the problem studied
here from the formidable task of modeling general protein-protein binding. In
docking, such constraints are challenging to incorporate into the scoring
functions and are approximated through separate steps to filter compatible PPI
poses and PROTAC geometries. While CGMD excludes many degrees of freedom from
the PROTAC, proteins, and solvent entropy, this effect of configurational
entropy on PPIs from PROTAC mediation can be directly captured.  Finally, we
calculate binding energies using alchemical methods, which circumvents the
computational challenge of directly sampling binding and unbinding events
between the PROTAC and proteins. We demonstrate the computational amenity of
an unconventional application of alchemical methods motivated by the PROTAC
systems, and take advantage of the physical interpretability of the CGMD +
alchemical approach to explore the principles of PROTAC binding cooperativity.

##  2 Methods

###  2.1 CGMD setup of PROTAC-protein complexes

The binary and ternary PROTAC-protein complexes are coarse-grained at two
resolutions to efficiently sample complex conformational changes while
retaining sufficient details for structural insight. Specifically, a major
focus of this work is to characterize the entropic effect of the length of
PROTACs on the strength of induced PPIs, necessitating modeling the PROTAC
linker at a higher resolution than the rest of the system. Proteins are
coarse-grained by mapping every three amino acids onto a large bead of
$\sigma=0.8$ nm diameter, which is approximately the Kuhn length of
polypeptides 43, 44, 45, 46. Binding moieties at the two ends of a PROTAC are
each represented by a large bead, whereas the linker region is modeled as a
Gaussian chain at the resolution of a PEG unit ($\sigma_{s}=0.35$ nm 47) or
three heavy atoms. Several experimental works that used flexible linear
linkers motivate our modeling approach for the PROTAC linker, including Chan
et al. 28 where an alkane linker was varied in step sizes of our linker beads
and Zorba et al. 29 where a PEG linker is modified at smaller length steps
such that linker lengths ranging from 1 to 6 $\sigma_{s}$ in our modeling
correspond to the PROTAC (1), (3), (5), (6), (8), and (10).

A minimal forcefield is used to describe the internal and interactive forces,
and a full description can be found in the Supporting Information (Section
S1). The three-dimensional structure of a protein is maintained by a bottom-up
fitted elastic network model (Fig. S2), which allows conformational
flexibility 48, 49. Protein beads can have additional properties to describe
PPIs beyond volume exclusion (Fig. S1). When modeling electrostatic
interactions, for example, a protein bead has the net charges of the triplet
of residues that it is coarse-grained from. PROTACs are modeled as Gaussian
polymers with volume exclusion, and the warhead beads are attached to the
binding pockets of proteins through harmonic springs. Modeling PROTAC
interactions beyond warhead binding is out of the scope of this work. Thus,
under current setup, PROTAC beads have 0 charge and no affinity to any other
beads.

The orientation between the E3 ligase and the target protein is initialized
such that the two binding pockets face each other, with a fully extended
PROTAC tethering in between (Fig. 1a). The binding moiety beads of PROTAC are
placed at the center of each binding pocket, which is defined by the residues
within 4 or 5 Å from the PROTAC warhead in experimental structures. Thus,
setting up the initial coordinates of a ternary complex requires the following
inputs: structures of each protein, residues at the two PROTAC binding
pockets, and the length of the PROTAC linker. To calculate the difference in
PROTAC binding energies due to PPIs, simulations of binary target/E3-PROTAC
complexes are also needed. Binary complexes are prepared by removing a protein
from the initialized ternary complex.

Figure 1: Schematic of the simulation setup for PROTAC-mediated complexes. (a)
The target-PROTAC-E3 ternary complex is initialized with a fully extended
PROTAC as drawn. The proteins are coarse-grained at the resolution of three
amino acids per bead, approximately $0.8$ nm. PROTAC warhead beads are
represented by beads of the same size, whereas the linker is coarse-grained at
a higher resolution. (b) The thermodynamic cycle shows that the $\Delta\Delta
G$ is obtained by measuring the free energy difference between target-PROTAC
binding with and without the E3 ($\Delta G^{\text{binary}}_{TP}-\Delta
G^{\text{ternary}}_{TP}$) or between PROTAC-E3 binding with and without the
target ($\Delta G^{\text{binary}}_{EP}-\Delta G^{\text{ternary}}_{EP}$). Under
the alchemical setup, $\Delta\Delta G$ can be alternatively quantified by the
free energy difference between the red vertical processes, which represent
coupling the target ($\Delta G^{\text{binary}}_{\text{couple\\_}T}-\Delta
G^{\text{ternary}}_{\text{couple\\_}T}$ in (c)) or the E3 ($\Delta
G^{\text{binary}}_{\text{couple\\_}E}-\Delta
G^{\text{ternary}}_{\text{couple\\_}E}$ in (d)) to the PROTAC and the PROTAC
pre-bound to the other protein. In the initial states in (c) and (d), the
dotted lines represent the target or the E3 whose interactions with the rest
of the system are turned off except for the harmonic constraints (black lines)
to the PROTAC warhead.

###  2.2 Thermodynamic framework of alchemical perturbation

The binding cooperativity of a PROTAC is mathematically defined as
$\exp{\left(\frac{\Delta\Delta G}{RT}\right)}$, where $R$ is the gas constant,
$T$ here refers to the temperature in the context of an energetic scale and
refers to the target protein elsewhere, $\Delta\Delta G=\Delta
G^{\text{binary}}_{TP}-\Delta G^{\text{ternary}}_{TP}$, and $\Delta
G^{\text{ternary}}_{TP}$ and $\Delta G^{\text{binary}}_{TP}$ are the free
energies of the PROTAC ($P$) binding to the target protein ($T$) with and
without the presence of the E3 ligase ($E$). Because of the thermodynamic
cycle (Fig. 1b), the same $\Delta\Delta G$ can be obtained from $\Delta
G^{\text{binary}}_{EP}-\Delta G^{\text{ternary}}_{EP}$. Favorable PPIs
stabilize the ternary complex and facilitate PROTAC binding to both proteins.
Thus, they lower $\Delta G^{\text{ternary}}_{TP}$ and $\Delta
G^{\text{ternary}}_{EP}$, which leads to larger $\Delta\Delta G$ and more
positive cooperativity.

Alchemical free energy calculation methods exploit alternative thermodynamic
cycles to obtain $\Delta\Delta G$ without simulating binding and unbinding
processes. For simplicity, in this work, all $\Delta\Delta G$s are calculated
using the cycle in Fig. 1c, which we describe in detail here, but one should
arrive at the same result using the mirroring cycle in Fig. 1d. By the
definition of a thermodynamic cycle, we have $\Delta
G^{\text{binary}}_{EP^{\prime}}-\Delta G^{\text{ternary}}_{EP^{\prime}}=\Delta
G^{\text{binary}}_{\text{couple\\_}T}-\Delta
G^{\text{ternary}}_{\text{couple\\_}T}$, where $\Delta
G^{\text{binary}}_{\text{couple\\_}T}$ and $\Delta
G^{\text{ternary}}_{\text{couple\\_}T}$ represent the free energies of
coupling $T$ to $P$ and to the target-PROTAC bound complex $EP$. In the
initial states of both coupling processes (vertical processes in red in Fig.
1c), $T$ is bound to $P$ but is a dummy molecule at an ideal state.
Specifically, multiple harmonic springs connect the binding pocket beads in
$T$ to the warhead bead of $P$, and $T$ itself is an elastic network model
consisting of only harmonic springs. All other interactions between $T$ and
the rest of the system – whether $P$ or $EP$ – are turned off. Coupling $T$
simply means turning on these inter-molecular interactions, while the binding
pocket springs remain unperturbed.

Attaching a dummy $T$ instead of having $T$ dissociated results in a
systematic error in the horizontal free energies of $EP$ binding ($\Delta
G^{\text{binary}}_{EP^{\prime}}$ and $\Delta G^{\text{ternary}}_{EP^{\prime}}$
in Fig. 1c) such that the $\Delta\Delta G$ is unaffected. This is because the
attachment of dummy $T$ occurs via only one bead on $P$, except which there
are no other forcefield terms involving both physically present beads and
dummy beads. In the configurational partition function, energy terms
describing the geometries of the physically present part of the system can
therefore be separated from the term involving the dummy $T$ and the
attachment junction. The latter term is the same whether the physically
present part is $P$ or $EP$, such that the unphysical contribution from
attaching dummy $T$ cancels out in $\Delta\Delta G$.

### Free energy calculations

Alchemically creating a protein from a dummy state to full coupling involves
turning on the interaction potentials between the protein and the rest of the
system in the forcefield. The interactions are turned on in stages by
sequentially scaling each kind of interaction potential using a coupling
parameter $\lambda$. Intramolecular potentials (e.g. the elastic network model
of each protein) and intermolecular potentials not perturbed at the current
stage are unaffected by the $\lambda$ scaling. For the electrostatic
potential, the start state (no electrostatics) and the end state (full
electrostatics) correspond to $\lambda_{\text{elec}}=0$ and $1$ respectively.
Intermediate states are interpolated such that the potential is defined as
$U_{\lambda_{\text{elec}}}=(1-\lambda_{\text{elec}})U_{\text{no\\_elec}}+\lambda_{\text{elec}}U_{\text{elec}}=\lambda_{\text{elec}}U_{\text{elec}}$.
For numerical stability, the electrostatic potential is only perturbed in the
presence of volume exclusion 50, 51, which is modeled by Weeks-Chandler-
Andersen (WCA) potential. To turn on Lennard-Jones (LJ) or variants of LJ
potentials (e.g. WCA), a soft-core scaling 52 with $\lambda_{\text{LJ}}$ is
used for numerical stability:

$U_{\lambda_{\text{LJ}}}\left(r_{ij}\right)=4\epsilon\lambda_{\text{LJ}}\left(\frac{1}{\left(\alpha\left(1-\lambda_{\text{LJ}}\right)+\left(\frac{r_{ij}}{\sigma_{ij}}\right)^{6}\right)^{2}}-\frac{1}{\alpha\left(1-\lambda_{\text{LJ}}\right)+\left(\frac{r_{ij}}{\sigma_{ij}}\right)^{6}}\right),$

where $\alpha=0.5$, $r_{ij}$ is the distance between beads $i$ and $j$, and
$\sigma_{ij}$ is the sum of the radii of beads $i$ and $j$. The number of
intermediate states and the spacing of the coupling parameter values depend on
the difficulty to obtain converged free energy calculations. For the
electrostatic potential, a linear pathway where $\lambda_{\text{elec}}$ ranges
from 0 to 1 with a step size of 0.125 is a simple and effective approach. For
LJ and related potentials, because most of the free energy changes occur near
the start state of $\lambda_{\text{LJ}}=0$ (Fig. 2b,c), we introduce
intermediate states at
$\lambda_{\text{LJ}}=0.005,0.01,0.015,0.02,0.04,0.06,0.08,0.1,0.2,0.3,0.5,0.7,$
and $0.9$.

The $\Delta G$ of turning on each kind of interaction is calculated using
thermodynamic integration (TI) 53, Bennett acceptance ratio (BAR) method 54
and the multi-state BAR (MBAR) method 55. TI and BAR/MBAR are distinct
formulations for free energy calculations, and we verify that these methods
converge to similar values. The system in CGMD is evolved using overdamped
Langevin dynamics with a diffusion coefficient of 253
$\text{nm}^{2}\mathbin{/}$s and a timestep of 30 ns for stable integration. At
each state, at least 64 trajectories of 6 s long are generated to sample the
conformations of the complexes. After collecting the samples from
trajectories, post-processing involves calculating $\frac{\partial
U}{\partial\lambda}$ and $\Delta U_{ij}$ for all $i,j=1,2,...,K$ states as
inputs for TI, BAR, and MBAR.

##  3 Results and discussion

###  3.1 Alchemical perturbation of protein domains is feasible with CGMD

The binding cooperativity of PROTAC due to PPIs is a unique challenge that
calls for an unconventional application of alchemical free energy calculation
methods. Alchemical methods are mainly used to determine the binding energies
between small-molecule ligands and proteins, and typically no more than 10
heavy atoms are perturbed for efficient and accurate calculations. In protein-
protein binding, recent applications and development focus on quantifying the
relative free energy changes from small-scale perturbations such as mutations
of single residues 56, 57, 58, 59, 60. To our knowledge, the only case that
alchemically calculates PPIs in a three-body setting compares how analogs of
inhibitors change aberrant multimerization of the HIV-1 integrase 61. Their
proposed thermodynamic framework involves calculating the relative free energy
difference by perturbing small molecules that directly participate at a fixed
PPI interface. This framework is more readily extendable to molecular glues
that modulate PPIs in a similar way. PROTACs, however, due to a more modular
design, are typically larger linear molecules. The flexibility of the linker
is often nontrivial, such that the two proteins cannot be kept bound at a
fixed interface. This configurational entropic concern necessitates an
unusually large perturbation at the scale of a protein rather than a small
molecule to calculate the binding cooperativity, testing the computational
limit of alchemical methods.

To explore the feasibility of the CG alchemical approach, we calculate the
free energy of turning on the steric repulsions between a target protein and a
PROTAC-E3 complex ($\Delta G\textsuperscript{ternary(sterics)}$) in the
absence of other inter-molecular potentials. We choose Bruton’s tyrosine
kinase (BTK) as the target (only the kinase domain modeled), CRBN as the E3,
and the PROTAC (10) from 29, which are respectively modeled by 87, 124, and 8
beads in the CG model. Together they form the largest target-PROTAC-E3 complex
simulated in this work. We compare the calculations using different
percentages of the simulation data collected in the time-forward and time-
reversed directions. The calculated values of $\Delta
G\textsuperscript{ternary(sterics)}$ plateau starting around the midpoint of
the simulation time, indicating numerical convergence (Fig. 2a). The time-
forward and -reversed estimations are within 1 standard deviation (std) at the
midpoint, and the time-reversed estimations remain stable after the midpoint.
The observed behavior of the estimates over time suggests that unequilibrated
samples at the beginning of the trajectories have been removed, and the
remaining frames sample from similar distributions rather than distinct
metastable states with slow transition rates 51.

Three methods, TI, BAR, and MBAR are used to separately estimate the free
energies. The accuracy of all three methods depends on the number and the
spacing of alchemical states. BAR and MBAR reweight conformations sampled from
one state by their probability in another state to estimate the free energy
differences. Having similar probability distributions between states, i.e.
phase space overlap, is therefore critical to the estimation. Unlike BAR/MBAR,
TI estimates the free energies by numerically integrating
$\langle\frac{\partial U}{\partial\lambda}\rangle$, the ensemble average of
the derivative of the potential energy $U$ along the alchemical pathway
defined by $\lambda$. Depending on the curvature of $\langle\frac{\partial
U}{\partial\lambda}\rangle$, choices of intermediate states specified by
$\lambda$ and the integration scheme together introduce integration errors in
addition to the statistical errors in estimating the ensemble average per
state. We choose an alchemical pathway that involves 12 intermediate states in
addition to the start and end states, such that $\Delta
G\textsuperscript{ternary(sterics)}=\sum_{i=1}^{13}\Delta
G_{\lambda_{i},\lambda_{i+1}}$, where $\Delta G_{\lambda_{i},\lambda_{i+1}}$
is the free energy of changing the WCA potential between neighboring states
$\lambda_{i}$ and $\lambda_{i+1}$. With a total of 14 states unevenly spaced,
the phase space overlap between neighboring states is sufficient (Fig. S3) for
efficient reweighting-based estimations. For TI, the trapezoid rule of
numerical integration is used for its simplicity and robustness. Although the
quadrature errors result in a slight overestimation of $\Delta
G\textsuperscript{ternary(sterics)}$, the $\partial U/\partial\lambda$ curve
is sufficiently smooth such that TI and MBAR largely agree. In addition to the
global agreement on $\Delta G\textsuperscript{ternary(sterics)}$, TI, BAR, and
MBAR also locally agree with each other on all $\Delta
G_{\lambda_{i},\lambda_{i+1}}$ along the alchemical pathway (Fig. 2c). We
emphasize that TI and BAR/MBAR rely on distinct types of input data and
processing procedures, and their consistency even at the most granular level
of calculations further validate our CG alchemical approach.

Analyses of estimations over simulation time and using different free energy
calculation methods indicate that convergence of perturbing a protein can be
achieved within reasonable computation time, significantly pushing the
boundaries of applying alchemical methods. As parallelization can be done over
the alchemical states and over trajectories for each state, the time to run
one trajectory is the main limiting factor in the wall-clock computation time
of applying our method. Criteria to determine how long a trajectory should be
run are described in the Supporting Information (Section S2). For this work,
depending on the size of the system, 3 - 14 CPU hours per trajectory of
ternary complexes are sufficient.

Figure 2: Calculation of $\Delta G\textsuperscript{ternary(sterics)}$ by
alchemical perturbation of BTK in the ternary complex of BTK-PROTAC (10)-CRBN.
(a) TI and MBAR both reach apparent convergence in the time-forward and time-
reversed directions with no pathological signs. The grey band in each panel
represents the final estimation using $100\%$ data $\pm 0.1$ $kT$ as a
threshold for error tolerance, where $k$ is the Boltzmann constant. (b) TI
estimation is shown as the blue area under the curve of $\langle\partial
U/\partial\lambda\rangle$. (c) TI, BAR, and MBAR agree for all intermediate
$\Delta G$s between adjacent states. All error bars of computational results
here and in subsequent figures represent $\pm 1$ std. Color coding for TI,
BAR, and MBAR results are the same in subsequent figures unless otherwise
stated.

###  3.2 Minimal forcefield captures entropic effects in PROTAC-mediated PPIs

Encouraged by the proof-of-concept calculations above for $\Delta
G\textsuperscript{ternary}$, we also calculate $\Delta
G\textsuperscript{binary}$ and complete our calculations for the $\Delta\Delta
G$ of the thermodynamic cycle. We follow the sign convention of $\Delta\Delta
G$ such that a positive value represents positive cooperativity. The BTK-CRBN
system modeled here has been experimentally shown to lack large cooperativity,
and introducing PROTACs in Hydrogen/Deuterium Exchange experiments didn’t
reveal significant profile changes that would indicate the presence of stable
PPIs. As the starting point for our method development, we focus on this
system due to its apparent simplicity and the availability of experimental
characterization over a large range of PROTAC linker lengths. We characterize
$\Delta\Delta G$ changes over PROTAC lengths because this relies on capturing
the fundamental physics of the tertiary interactions rather than sequence- or
conformation-specific properties.

Two forcefield setups are used to describe PPIs and the resulted $\Delta\Delta
G$ trends over PROTAC linker lengths are compared. In the first setup, we
calculate the baseline $\Delta\Delta G$ in the absence of PPIs other than
volume exclusion. In the second setup, nonspecific attractions between BTK and
CRBN beads are added and explored at two strengths. The intrinsic PPIs without
PROTAC mediation should be weak such that in the limit of infinite linker
length the $\Delta\Delta G$ is negligible. The attenuation of weak PPIs with
increasing PROTAC linker lengths originates from configurational entropy. As
the PROTAC becomes longer, it experiences a greater loss of configurational
freedom upon binding to proteins to induce PPIs, incurring an entropic cost.
We examine this configurational entropic effect by modeling $\Delta\Delta G$
at linkers ranging from 1 to 6 beads ($\sigma_{s}$) long, which correspond to
approximately 3.5 Å to 21 Å.

In the first setup, the steric cores of the proteins should penalize PROTAC
binding and result in negative cooperativities. This is because some
conformations that are accessible to the PROTAC in a binary PROTAC-protein
complex become inaccessible in the ternary complex due to steric clashes. As
the linker length increases and steric clashes are attenuated, the
cooperativity should become less negative. We verify that such a monotonically
increasing trend of negative $\Delta\Delta G$ is obtained in our model (Fig.
3). Steric penalties on $\Delta\Delta G$ are most obvious at the region of
short linker lengths (1-3 beads), after which the benefit from extending the
linker length becomes increasingly marginal, and we expect that beyond the
simulated window of linker lengths, $\Delta\Delta G$ will eventually plateau
near 0. This $\Delta\Delta G$ trend is consistent with a recent effort to
tabulate PROTAC linker length structure-activity relationships (SAR), which
suggests that steric clashes at short linker lengths often result in a steep
decrease in activity 38.

After validating the baseline trend, we next examine how the cooperativity
trend is changed by the addition of favorable PPIs through LJ potentials.
Increasing the well depth of LJ ($\epsilon_{\text{LJ}}$) increases the
strength of this nonspecific attraction, which is kept weak (Fig. S1) to
approximate van der Waals forces. At the attraction strength of
$\epsilon_{\text{LJ}}=0.125$ $kT$, the $\Delta\Delta G$ curve is elevated
compared to the previous curve without attraction (Fig. 3a), as favorable PPIs
are expected to enhance cooperativity. Nevertheless, at this attraction
strength, steric penalties still dominate and $\Delta\Delta G$s remains
negative. Even though adding an LJ potential brings an additional penalty when
beads overlap, shorter PROTACs still benefit more from the attractive part of
LJ than longer PROTACs, resulting in a flatter $\Delta\Delta G$ trend as
compared with the purely repulsive interactions.

Figure 3: The $\Delta\Delta G$ trends over PROTAC linker lengths change with
varying strengths of non-specific attraction between proteins. All
calculations shown are obtained using MBAR, and results using TI and BAR are
superimposed in Figure S5. The $\Delta\Delta G$ trends are calculated for two
target-E3 pairs, BTK-CRBN (a) and BRD4BD2-VHL (b). Solid lines represent the
baseline $\Delta\Delta G$ trends where only volume exclusion is modeled
between the two proteins, and dotted lines show the trends where nonspecific
attractions are added. The strengths ($\epsilon_{\text{LJ}}$) of attractions
are indicated by different colors. Higher $\epsilon_{\text{LJ}}$ represents
stronger attractions, and the baselines can also be considered as results at
$\epsilon_{\text{LJ}}=0$. Results at $\epsilon_{\text{LJ}}=0.125$ and $0.2$
$kT$ are plotted for BTK-CRBN and results at $\epsilon_{\text{LJ}}=$0.125,
0.15, 0.175, 0.2, and 0.25 $kT$ are plotted for BRD4BD2-VHL.

An appropriate combination of repulsive and attractive forces may generate a
non-monotonic $\Delta\Delta G$ trend, such that intermediate linker lengths
promote optimal cooperativity by minimizing steric clashes while maximally
sampling attractive PPIs38. As the attraction strength increases to
$\epsilon_{\text{LJ}}=0.2$ $kT$, intermediate-length PROTACs exhibit not only
positive $\Delta\Delta G$s but the values can be comparable and even slightly
higher than that of the longest PROTAC (Fig. 3a). Within the limited window of
linker lengths, only the initial part of the decaying tail of a non-monotonic
$\Delta\Delta G$ trend is observed. We expect that beyond the simulated window
of linker lengths, configurational entropic penalties will continue driving
$\Delta\Delta G$ down towards 0.

Experimentally, the linker length at 3 beads uniquely enables weak positive
cooperativity for BTK-CRBN, whereas our results at $\epsilon_{\text{LJ}}=0.2$
$kT$ remain biased towards favoring longer linkers and are not as sensitive to
linker length changes. To see whether these characteristics are specific to
the choice of the system, we then examine the $\Delta\Delta G$ trends for a
different system (Fig. 3b), BRD4BD2-VHL, where experimentally, the linker
length at 3 beads can also optimize the cooperativity 28. Due to the smaller
size of the system, we can afford to calculate $\Delta\Delta G$s at three more
attraction strengths. Similar to BTK-CRBN, in the absence of attractions,
negative $\Delta\Delta G$ monotonically increases over the linker length, and
adding nonspecific attractions results in flatter and higher $\Delta\Delta G$
curves. Within the narrow window of short linker lengths, scanning the
attractive strength $\epsilon_{\text{LJ}}$ from 0.125 to 0.25 $kT$, however,
does not recapitulate the optimal linker length at 3 beads. This result
suggests that enhancing nonspecific attractions in the minimal model is
insufficient to compensate for the steric penalties while remaining sensitive
to entropic penalties from the linker length.

We demonstrate that the minimal CG model directly captures configurational
entropic effects on weak nonspecific PPIs through analyzing $\Delta\Delta G$
trends over PROTAC linker lengths. Beyond this entropic effect, combining
repulsive and attractive interactions at various strengths changes the
behaviors of cooperativity trends and can shift the optimal linker length, as
shown in BTK-CRBN. Nevertheless, chemically specific interactions or specific
sampling of certain PPIs is needed to model optimal positive cooperativity at
an experimentally relevant range and resolution of PROTAC linker lengths.

### Electrostatics in PROTAC-mediated PPIs exhibit plasticity

As a step towards more realistic modeling of cooperativity, we seek chemically
specific PPIs to include and further explore the BRD4BD2-VHL system due to the
availability of experimental structural information. Crystal structures of the
ternary complexes have revealed specific interactions that are proposed as the
molecular basis for the observed positive cooperativity and selectivity
against other structural homologs.8, 62 As shown in the previous subsection,
these interactions between proteins cannot be approximated by nonspecific
attractions that contribute to the cooperativity with low sensitivity to
linker length and no protein sequence dependence.

The structural findings such as salt bridges at the PPI interface and the
mutational studies involving charged residues on BRD4BD2 and homologs8
motivate us to approach chemical specificity through modeling electrostatic
interactions. As CGMD uses an implicit solvent, we choose the Debye-Hückel
(DH) potential to describe electrostatics in consideration of screening
effects under physiological conditions. Within the BRD4BD2-VHL system,
incorporating charges of protein beads results in a monotonic trend of
negative $\Delta\Delta G$s with increasing linker length, (Fig. 4a) similar to
the baseline obtained using steric repulsions only (Fig. 3b). Since charges
are perturbed separately in $\Delta\Delta G$ calculations for numeric
stability, in the following discussions, we further investigate our
$\Delta\Delta G$ results by isolating the final stage ($\Delta
G\textsuperscript{ternary(charges)}$) in which charges are turned on in the
presence of sterics.

Breaking down the $\Delta\Delta G$s by each energy component shows that at all
three linker lengths, $\Delta G\textsuperscript{ternary(charges)}$ is slightly
negative, indicating a mildly favorable process, but the penalty from steric
repulsions overwhelmingly dominates electrostatic contributions by an order of
magnitude (Fig. 4c). As PROTAC linker length increases from 2 to 4 beads, the
contribution from $\Delta G\textsuperscript{ternary(charges)}$ monotonically
diminishes. We consider the possibility that the screening of charges is too
strong to model more favorable PPIs and tune the screening parameter in the DH
potential at the linker length of 3 beads. However, because both the target
protein and the E3 ligase have net positive charges, significantly weakening
the screening strength leads to a much more unfavorable $\Delta
G\textsuperscript{ternary(charges)}$ (Fig. 4c). It is also possible that our
level of coarse-graining loses the spatial resolution required for this system
to capture detailed interactions like salt bridge formation as observed in the
crystal structures8, 62.

In addition to the small contribution to $\Delta\Delta G$, $\Delta
G\textsuperscript{ternary(charges)}$ itself exhibits plasticity because
conformational sampling at the stage of charge perturbation in alchemical free
energy calculations is biased by the potentials turned on in previous stages.
The presence of steric repulsions combined with nonspecific attractions at the
strength of $\epsilon_{\text{LJ}}=0.2$ $kT$, for example, has doubled the
$\Delta G\textsuperscript{ternary(charges)}$ obtained at the linker length of
3 beads without nonspecific attractions (Fig. 4c). Interestingly, this change
in $\Delta G\textsuperscript{ternary(charges)}$ is on top of the favorable
contribution from nonspecific attractions in the previous calculation stage
($\Delta G\textsuperscript{ternary(other)}$) before the inclusion of protein
charges. For this particular ternary complex, nonspecific attractions and
electrostatic interactions work synergistically.

Our dissection of the electrostatic component in $\Delta\Delta G$ under
different simulation setups suggests that a more holistic parameterization is
needed to accurately evaluate chemically specific PPIs. For BRD4BD2-VHL,
incorporating hydrophobic interactions will be of particular interest as there
is stacking of hydrophobic residues at the PPI interface observed in the
crystal structures.8, 62 Hydrophobic interactions may also introduce non-
additive free-energy contributions with electrostatics in a similar manner
seen with the nonspecific attractions. It is also worth noting that the
favorable PPIs revealed by crystal structures are enabled by PROTACs using a
JQ1 warhead, which imposes a different linker attachment angle (i.e. exit
vector) from an I-BET726 warhead (Fig. S7).28 Our current forcefield does not
model the PROTAC linker with angular terms to specify the exit vectors, which
leads to a $\Delta\Delta G$ trend that matches well with the worse-performing
I-BET726 set of PROTACs (Fig. 4a). As rigidifying PROTACs is a common strategy
to optimize the cooperativity by entropically enhancing certain PPIs62, 30,
parameterizing linker conformations will improve modeling the specificity in
PROTAC-mediated PPIs.

Figure 4:  Electrostatic contributions to the cooperativity in the BRD4BD2-VHL
system are small and context-dependent. All calculations shown are obtained
using MBAR, and results using TI and BAR are shown in Figure S6. (a)
Calculations of $\Delta\Delta G$s over PROTAC linker lengths are shown with
the experimental measurements 28 (black) converted to our units. Experimental
results at 2, 3, and 4 linker beads correspond to MZ4, MZ1, and MZ2 for
PROTACs using JQ1 warhead and MZP-61, MZP-54, MZP-55 for PROTACs using
I-BET726 warhead. (b) Waterfall plot breakdown of $\Delta\Delta G$
calculations. At each linker length, bars in each triplet correspond to
$\Delta G\textsuperscript{binary}$ (grey), $-\Delta
G\textsuperscript{ternary(other)}$ (light purple), and $-\Delta
G\textsuperscript{ternary(charges)}$ (turquoise), and are arranged in a
cumulative manner such that the end position marks the resulted $\Delta\Delta
G$ (orange). $\Delta G\textsuperscript{ternary(other)}$ denotes the free
energy change of turning on interaction energy components other than the
electrostatics, which only include steric repulsions in this panel. (c)
$\Delta\Delta G$ breakdowns at linker length 3 under different forcefield
parameterizations are superimposed for comparison. Reducing the screening
effect by ten-fold (charges*) significantly increases $\Delta
G\textsuperscript{ternary(charges)}$ (cyan), which leads to a very negative
$\Delta\Delta G$. Introducing non-specific attractions
($\epsilon_{\text{LJ}}=0.2$ $kT$) not only reduces $\Delta
G\textsuperscript{ternary(other)}$ (dark purple) but also doubles $\Delta
G\textsuperscript{ternary(charges)}$ (steel blue), resulting in a positive
$\Delta\Delta G$.

##  4 Conclusions

We explore a novel computational approach to model the binding cooperativity
of PROTACs by combining CGMD and alchemical free energy calculations. The
plasticity of PROTAC-mediated PPIs motivates an unconventional application of
alchemical methods at a perturbation scale that is rarely attempted. We show
that with coarse-graining, converged estimates from various free energy
calculation methods are attainable within a reasonable amount of computation
time. Our results expand the possibility of more creative use of alchemical
methods. The feasibility and efficiency of the CG alchemical approach enable
us to probe multiple energy components under the alchemical framework and
characterize how PROTAC linker lengths modulate PPIs under different setups to
produce distinct cooperativity trends. In addition to validating the benefit
of using long linkers to avoid steric clashes, we demonstrate with a simple
addition of nonspecific attractions between BTK and CRBN that the binding
cooperativity can be promoted by shortening the PROTAC linker. Our minimal
model is capable of unveiling such changes in cooperativity that are rooted in
the configurational freedom of the ternary complexes rather than chemical
properties.

Quantitative modeling of the cooperativity, however, remains difficult due to
the lack of specificity in the minimal model. Previous studies have recognized
the challenges brought by non-native PROTAC-mediated PPIs that are often weak,
transient, and pliable, and have called for a paradigm shift towards an
ensemble-based characterization beyond a handful of docked or crystal poses.
9, 34, 39. While thermodynamic properties such as the binding cooperativity
are inherently ensemble-based, we note that both accurate sampling of PPI
conformations according to chemical properties and efficient computation to
sample a diverse set of conformations are important for calculations.
Currently, tuning the strength of nonspecific attractions cannot approximate
favorable PPIs while retaining sensitivity to entropic constraints from the
PROTAC linker length. Simply adding electrostatic interactions based on amino
acid charges proved insufficient to capture the cooperativity trend enabled by
JQ1-based PROTACs in BRD4BD2-VHL. Additional parameterizations are needed to
capture chemically specific PPIs.

Two main avenues are worth exploring for future improvement of our method –
PROTAC linker conformations and protein sequence-dependence. Among a myriad of
PROTAC properties 6 that we leave out, structural features such as the exit
vector 28 and the linker rigidity 62, 30 in addition to the linker length can
both entropically constrain the sampling of PPIs. Meanwhile, energy components
of PPIs other than electrostatic interactions, notably the hydrophobic
effects, are currently overlooked. Different energy components may have non-
additive effects in optimizing the absolute cooperativity and relative
cooperativities between target homologs such as BRD4BD2 and BRD4BD1. Although
coarse-graining enables efficient computation, parameterization for both
directions of forcefield development will be a major hurdle to overcome. This
can be bottom-up using shorter-timescale higher-resolution simulations,
similar to that of the CG ENM (Fig. S2) in this work. A top-down fitting might
also become possible with rapidly growing experimental studies that develop
platforms 63 for empirical SAR of PROTAC linkerology 64, 65 or leverage
promiscuous PROTACs and target homologs and mutants to investigate the
molecular basis of specificity 66.

Description of the forcefield terms in CGMD, parameterization of CG ENM,
analysis of phase space overlap in alchemical free energy calculations of the
BTK-PROTAC (10)-CRBN complex, description of post-processing equilibrated and
statistically de-correlated samples from CGMD trajectories for free energy
calculations, convergence among TI, BAR, and MBAR for the results shown in
Figure 3 and 4, and crystal structures showing the exit vector difference
between JQ1 and I-BET726.

H.M. thanks William M. Clemons, Jr., Daniel Jacobson, Tomislav Begušić,
Xuecheng Tao, Marta Gonzalvo, and Lixue Cheng for comments on the manuscript,
and Zhen-Gang Wang and Christopher J. Balzer for technical discussions. We
gratefully acknowledge support from the National Institutes of Health (NIH)
R01GM138845 (8877_CIT, subaward), Amgen Chem-Bio-Engineering Award (CBEA), and
DeLogi Trust Science and Technology Grant. This work used the Extreme Science
and Engineering Discovery Environment (XSEDE) Bridges computer at the
Pittsburgh Supercomputing Center through allocation MCB160013 67. XSEDE is
supported by National Science Foundation grant number ACI-1548562. This work
also used computational resources from the Resnick High Performance Computing
Center, a facility supported by Resnick Sustainability Institute at the
California Institute of Technology.

## References

  * An and Fu 2018 An, S.; Fu, L. Small-molecule PROTACs: An emerging and promising approach for the development of targeted therapy drugs. _EBioMedicine_ 2018, _36_ , 553–562
  * Burslem and Crews 2020 Burslem, G. M.; Crews, C. M. Proteolysis-Targeting Chimeras as Therapeutics and Tools for Biological Discovery. _Cell_ 2020, _181_ , 102–114
  * Sakamoto et al. 2001 Sakamoto, K. M.; Kim, K. B.; Kumagai, A.; Mercurio, F.; Crews, C. M.; Deshaies, R. J. Protacs: Chimeric molecules that target proteins to the Skp1–Cullin–F box complex for ubiquitination and degradation. _Proceedings of the National Academy of Sciences_ 2001, _98_ , 8554–8559, Publisher: National Academy of Sciences Section: Biological Sciences
  * Mullard 2021 Mullard, A. Targeted protein degraders crowd into the clinic. _Nature Reviews Drug Discovery_ 2021, _20_ , 247–250, Bandiera_abtest: a Cg_type: News Number: 4 Publisher: Nature Publishing Group
  * Maniaci and Ciulli 2019 Maniaci, C.; Ciulli, A. Bifunctional chemical probes inducing protein–protein interactions. _Current Opinion in Chemical Biology_ 2019, _52_ , 145–156
  * Troup et al. 2020 Troup, R. I.; Fallan, C.; Baud, M. G. J. Current strategies for the design of PROTAC linkers: a critical review. _Exploration of Targeted Anti-tumor Therapy_ 2020, _1_ , 273–312, Publisher: Open Exploration
  * Alabi and Crews 2021 Alabi, S.; Crews, C. Major Advances in Targeted Protein Degradation: PROTACs, LYTACs, and MADTACs. _Journal of Biological Chemistry_ 2021, 100647
  * Gadd et al. 2017 Gadd, M. S.; Testa, A.; Lucas, X.; Chan, K.-H.; Chen, W.; Lamont, D. J.; Zengerle, M.; Ciulli, A. Structural basis of PROTAC cooperative recognition for selective protein degradation. _Nature Chemical Biology_ 2017, _13_ , 514–521, Number: 5 Publisher: Nature Publishing Group
  * Nowak et al. 2018 Nowak, R. P.; DeAngelo, S. L.; Buckley, D.; He, Z.; Donovan, K. A.; An, J.; Safaee, N.; Jedrychowski, M. P.; Ponthier, C. M.; Ishoey, M. et al. Plasticity in binding confers selectivity in ligand-induced protein degradation. _Nature Chemical Biology_ 2018, _14_ , 706–714, Number: 7 Publisher: Nature Publishing Group
  * Smith et al. 2019 Smith, B. E.; Wang, S. L.; Jaime-Figueroa, S.; Harbin, A.; Wang, J.; Hamman, B. D.; Crews, C. M. Differential PROTAC substrate specificity dictated by orientation of recruited E3 ligase. _Nature Communications_ 2019, _10_ , 131, Number: 1 Publisher: Nature Publishing Group
  * Riching et al. 2018 Riching, K. M.; Mahan, S.; Corona, C. R.; McDougall, M.; Vasta, J. D.; Robers, M. B.; Urh, M.; Daniels, D. L. Quantitative Live-Cell Kinetic Degradation and Mechanistic Profiling of PROTAC Mode of Action. _ACS Chemical Biology_ 2018, _13_ , 2758–2770, Publisher: American Chemical Society
  * Roy et al. 2019 Roy, M. J.; Winkler, S.; Hughes, S. J.; Whitworth, C.; Galant, M.; Farnaby, W.; Rumpel, K.; Ciulli, A. SPR-Measured Dissociation Kinetics of PROTAC Ternary Complexes Influence Target Degradation Rate. _ACS Chemical Biology_ 2019, _14_ , 361–368, Publisher: American Chemical Society
  * Farnaby et al. 2019 Farnaby, W.; Koegl, M.; Roy, M. J.; Whitworth, C.; Diers, E.; Trainor, N.; Zollman, D.; Steurer, S.; Karolyi-Oezguer, J.; Riedmueller, C. et al. BAF complex vulnerabilities in cancer demonstrated via structure-based PROTAC design. _Nature Chemical Biology_ 2019, _15_ , 672–680, Number: 7 Publisher: Nature Publishing Group
  * Du et al. 2019 Du, X.; Volkov, O. A.; Czerwinski, R. M.; Tan, H.; Huerta, C.; Morton, E. R.; Rizzi, J. P.; Wehn, P. M.; Xu, R.; Nijhawan, D. et al. Structural Basis and Kinetic Pathway of RBM39 Recruitment to DCAF15 by a Sulfonamide Molecular Glue E7820. _Structure_ 2019, _27_ , 1625–1633.e3
  * Lai et al. 2016 Lai, A. C.; Toure, M.; Hellerschmied, D.; Salami, J.; Jaime-Figueroa, S.; Ko, E.; Hines, J.; Crews, C. M. Modular PROTAC Design for the Degradation of Oncogenic BCR-ABL. _Angewandte Chemie International Edition_ 2016, _55_ , 807–810, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.201507634
  * Bondeson et al. 2018 Bondeson, D. P.; Smith, B. E.; Burslem, G. M.; Buhimschi, A. D.; Hines, J.; Jaime-Figueroa, S.; Wang, J.; Hamman, B. D.; Ishchenko, A.; Crews, C. M. Lessons in PROTAC Design from Selective Degradation with a Promiscuous Warhead. _Cell Chemical Biology_ 2018, _25_ , 78–87.e5
  * Donovan et al. 2020 Donovan, K. A.; Ferguson, F. M.; Bushman, J. W.; Eleuteri, N. A.; Bhunia, D.; Ryu, S.; Tan, L.; Shi, K.; Yue, H.; Liu, X. et al. Mapping the Degradable Kinome Provides a Resource for Expedited Degrader Development. _Cell_ 2020, _183_ , 1714–1731.e10
  * Spradlin et al. 2019 Spradlin, J. N.; Hu, X.; Ward, C. C.; Brittain, S. M.; Jones, M. D.; Ou, L.; To, M.; Proudfoot, A.; Ornelas, E.; Woldegiorgis, M. et al. Harnessing the anti-cancer natural product nimbolide for targeted protein degradation. _Nature Chemical Biology_ 2019, _15_ , 747–755, Number: 7 Publisher: Nature Publishing Group
  * Ward et al. 2019 Ward, C. C.; Kleinman, J. I.; Brittain, S. M.; Lee, P. S.; Chung, C. Y. S.; Kim, K.; Petri, Y.; Thomas, J. R.; Tallarico, J. A.; McKenna, J. M. et al. Covalent Ligand Screening Uncovers a RNF4 E3 Ligase Recruiter for Targeted Protein Degradation Applications. _ACS Chemical Biology_ 2019, _14_ , 2430–2440, Publisher: American Chemical Society
  * Zhang et al. 2019 Zhang, X.; Crowley, V. M.; Wucherpfennig, T. G.; Dix, M. M.; Cravatt, B. F. Electrophilic PROTACs that degrade nuclear proteins by engaging DCAF16. _Nature Chemical Biology_ 2019, _15_ , 737–746, Number: 7 Publisher: Nature Publishing Group
  * Kuljanin et al. 2021 Kuljanin, M.; Mitchell, D. C.; Schweppe, D. K.; Gikandi, A. S.; Nusinow, D. P.; Bulloch, N. J.; Vinogradova, E. V.; Wilson, D. L.; Kool, E. T.; Mancias, J. D. et al. Reimagining high-throughput profiling of reactive cysteines for cell-based screening of large electrophile libraries. _Nature Biotechnology_ 2021, _39_ , 630–641, Number: 5 Publisher: Nature Publishing Group
  * Li et al. 2008 Li, W.; Bengtson, M. H.; Ulbrich, A.; Matsuda, A.; Reddy, V. A.; Orth, A.; Chanda, S. K.; Batalov, S.; Joazeiro, C. A. P. Genome-Wide and Functional Annotation of Human E3 Ubiquitin Ligases Identifies MULAN, a Mitochondrial E3 that Regulates the Organelle’s Dynamics and Signaling. _PLOS ONE_ 2008, _3_ , e1487, Publisher: Public Library of Science
  * Jevtić et al. 2021 Jevtić, P.; Haakonsen, D. L.; Rapé, M. An E3 ligase guide to the galaxy of small-molecule-induced protein degradation. _Cell Chemical Biology_ 2021,
  * Scholes et al. 2021 Scholes, N. S.; Mayor-Ruiz, C.; Winter, G. E. Identification and selectivity profiling of small-molecule degraders via multi-omics approaches. _Cell Chemical Biology_ 2021,
  * Huang et al. 2018 Huang, H.-T.; Dobrovolsky, D.; Paulk, J.; Yang, G.; Weisberg, E. L.; Doctor, Z. M.; Buckley, D. L.; Cho, J.-H.; Ko, E.; Jang, J. et al. A Chemoproteomic Approach to Query the Degradable Kinome Using a Multi-kinase Degrader. _Cell Chemical Biology_ 2018, _25_ , 88–99.e6
  * Rodriguez-Rivera and Levi 2021 Rodriguez-Rivera, F. P.; Levi, S. M. Unifying Catalysis Framework to Dissect Proteasomal Degradation Paradigms. _ACS Central Science_ 2021, Publisher: American Chemical Society
  * Maniaci et al. 2017 Maniaci, C.; Hughes, S. J.; Testa, A.; Chen, W.; Lamont, D. J.; Rocha, S.; Alessi, D. R.; Romeo, R.; Ciulli, A. Homo-PROTACs: bivalent small-molecule dimerizers of the VHL E3 ubiquitin ligase to induce self-degradation. _Nature Communications_ 2017, _8_ , 830, Number: 1 Publisher: Nature Publishing Group
  * Chan et al. 2018 Chan, K.-H.; Zengerle, M.; Testa, A.; Ciulli, A. Impact of Target Warhead and Linkage Vector on Inducing Protein Degradation: Comparison of Bromodomain and Extra-Terminal (BET) Degraders Derived from Triazolodiazepine (JQ1) and Tetrahydroquinoline (I-BET726) BET Inhibitor Scaffolds. _Journal of Medicinal Chemistry_ 2018, _61_ , 504–513, Publisher: American Chemical Society
  * Zorba et al. 2018 Zorba, A.; Nguyen, C.; Xu, Y.; Starr, J.; Borzilleri, K.; Smith, J.; Zhu, H.; Farley, K. A.; Ding, W.; Schiemer, J. et al. Delineating the role of cooperativity in the design of potent PROTACs for BTK. _Proceedings of the National Academy of Sciences_ 2018, _115_ , E7285–E7292, ISBN: 9781803662114 Publisher: National Academy of Sciences Section: PNAS Plus
  * Schiemer et al. 2021 Schiemer, J.; Horst, R.; Meng, Y.; Montgomery, J. I.; Xu, Y.; Feng, X.; Borzilleri, K.; Uccello, D. P.; Leverett, C.; Brown, S. et al. Snapshots and ensembles of BTK and cIAP1 protein degrader ternary complexes. _Nature Chemical Biology_ 2021, _17_ , 152–160, Number: 2 Publisher: Nature Publishing Group
  * Drummond et al. 2020 Drummond, M. L.; Henry, A.; Li, H.; Williams, C. I. Improved Accuracy for Modeling PROTAC-Mediated Ternary Complex Formation and Targeted Protein Degradation via New In Silico Methodologies. _Journal of Chemical Information and Modeling_ 2020, _60_ , 5234–5254, Publisher: American Chemical Society
  * Zaidman et al. 2020 Zaidman, D.; Prilusky, J.; London, N. PRosettaC: Rosetta Based Modeling of PROTAC Mediated Ternary Complexes. _Journal of Chemical Information and Modeling_ 2020, _60_ , 4894–4903, Publisher: American Chemical Society
  * Weng et al. 2021 Weng, G.; Li, D.; Kang, Y.; Hou, T. Integrative Modeling of PROTAC-Mediated Ternary Complexes. _Journal of Medicinal Chemistry_ 2021, _64_ , 16271–16281, Publisher: American Chemical Society
  * Bai et al. 2021 Bai, N.; Miller, S. A.; Andrianov, G. V.; Yates, M.; Kirubakaran, P.; Karanicolas, J. Rationalizing PROTAC-Mediated Ternary Complex Formation Using Rosetta. _Journal of Chemical Information and Modeling_ 2021, _61_ , 1368–1382, Publisher: American Chemical Society
  * Bai et al. 2022 Bai, N.; Riching, K. M.; Makaju, A.; Wu, H.; Acker, T. M.; Ou, S.-C.; Zhang, Y.; Shen, X.; Bulloch, D.; Rui, H. et al. Modeling the CRL4A ligase complex to predict target protein ubiquitination induced by cereblon-recruiting PROTACs. _Journal of Biological Chemistry_ 2022, 101653
  * Moreira et al. 2010 Moreira, I. S.; Fernandes, P. A.; Ramos, M. J. Protein–protein docking dealing with the unknown. _Journal of Computational Chemistry_ 2010, _31_ , 317–342, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.21276
  * Gromiha et al. 2017 Gromiha, M. M.; Yugandhar, K.; Jemimah, S. Protein–protein interactions: scoring schemes and binding affinity. _Current Opinion in Structural Biology_ 2017, _44_ , 31–38
  * Bemis et al. 2021 Bemis, T. A.; La Clair, J. J.; Burkart, M. D. Unraveling the Role of Linker Design in Proteolysis Targeting Chimeras. _Journal of Medicinal Chemistry_ 2021, Publisher: American Chemical Society
  * Eron et al. 2021 Eron, S. J.; Huang, H.; Agafonov, R. V.; Fitzgerald, M. E.; Patel, J.; Michael, R. E.; Lee, T. D.; Hart, A. A.; Shaulsky, J.; Nasveschuk, C. G. et al. Structural Characterization of Degrader-Induced Ternary Complexes Using Hydrogen–Deuterium Exchange Mass Spectrometry and Computational Modeling: Implications for Structure-Based Design. _ACS Chemical Biology_ 2021, Publisher: American Chemical Society
  * Dixon et al. 2021 Dixon, T.; MacPherson, D.; Mostofian, B.; Dauzhenka, T.; Lotz, S.; McGee, D.; Shechter, S.; Shrestha, U. R.; Wiewiora, R.; McDargh, Z. A. et al. _Atomic-Resolution Prediction of Degrader-mediated Ternary Complex Structures by Combining Molecular Simulations with Hydrogen Deuterium Exchange_ ; 2021; p 2021.09.26.461830, Section: New Results Type: article
  * Li et al. 2022 Li, W.; Zhang, J.; Guo, L.; Wang, Q. Importance of Three-Body Problems and Protein–Protein Interactions in Proteolysis-Targeting Chimera Modeling: Insights from Molecular Dynamics Simulations. _Journal of Chemical Information and Modeling_ 2022, Publisher: American Chemical Society
  * Liao et al. 2022 Liao, J.; Nie, X.; Unarta, I. C.; Ericksen, S. S.; Tang, W. In Silico Modeling and Scoring of PROTAC-Mediated Ternary Complex Poses. _Journal of Medicinal Chemistry_ 2022, _65_ , 6116–6132, Publisher: American Chemical Society
  * Niesen et al. 2017 Niesen, M. J. M.; Wang, C. Y.; Lehn, R. C. V.; Miller, T. F. Structurally detailed coarse-grained model for Sec-facilitated co-translational protein translocation and membrane integration. _PLOS Computational Biology_ 2017, _13_ , e1005427, Publisher: Public Library of Science
  * Zhang and Miller 2012 Zhang, B.; Miller, T. F. Long-Timescale Dynamics and Regulation of Sec-Facilitated Protein Translocation. _Cell Reports_ 2012, _2_ , 927–937, Publisher: Elsevier
  * Hanke et al. 2010 Hanke, F.; Serr, A.; Kreuzer, H. J.; Netz, R. R. Stretching single polypeptides: The effect of rotational constraints in the backbone. _EPL (Europhysics Letters)_ 2010, _92_ , 53001, Publisher: IOP Publishing
  * Staple et al. 2008 Staple, D. B.; Payne, S. H.; Reddin, A. L. C.; Kreuzer, H. J. Model for Stretching and Unfolding the Giant Multidomain Muscle Protein Using Single-Molecule Force Spectroscopy. _Physical Review Letters_ 2008, _101_ , 248301, Publisher: American Physical Society
  * Charmainne Cruje and Devika B Chithrani 2014 Charmainne Cruje,; Devika B Chithrani, Polyethylene Glycol Density and Length Affects Nanoparticle Uptake by Cancer Cells. _Journal of Nanomedicine Research_ 2014, _Volume 1_ , Publisher: MedCrave Publishing
  * Lezon et al. 2009 Lezon, T. R.; Shrivastava, I. H.; Yang, Z.; Bahar, I. _Handbook on Biological Networks_ ; World Scientific Lecture Notes in Complex Systems Volume 10; WORLD SCIENTIFIC, 2009; Vol. Volume 10; pp 129–158
  * Ricardo Batista et al. 2010 Ricardo Batista, P.; Herbert Robert, C.; Maréchal, J.-D.; Ben Hamida-Rebaï, M.; Geraldo Pascutti, P.; Mascarello Bisch, P.; Perahia, D. Consensus modes, a robust description of protein collective motions from multiple-minima normal mode analysis—application to the HIV-1 protease. _Physical Chemistry Chemical Physics_ 2010, _12_ , 2850–2859, Publisher: Royal Society of Chemistry
  * Pohorille et al. 2010 Pohorille, A.; Jarzynski, C.; Chipot, C. Good Practices in Free-Energy Calculations. _The Journal of Physical Chemistry B_ 2010, _114_ , 10235–10253, Publisher: American Chemical Society
  * Klimovich et al. 2015 Klimovich, P. V.; Shirts, M. R.; Mobley, D. L. Guidelines for the analysis of free energy calculations. _Journal of Computer-Aided Molecular Design_ 2015, _29_ , 397–411
  * Mey et al. 2020 Mey, A. S. J. S.; Allen, B.; Macdonald, H. E. B.; Chodera, J. D.; Kuhn, M.; Michel, J.; Mobley, D. L.; Naden, L. N.; Prasad, S.; Rizzi, A. et al. Best Practices for Alchemical Free Energy Calculations. _Living Journal of Computational Molecular Science_ 2020, _2_ , arXiv: 2008.03067
  * Kirkwood 1935 Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. _The Journal of Chemical Physics_ 1935, _3_ , 300–313, Publisher: American Institute of Physics
  * Bennett 1976 Bennett, C. H. Efficient estimation of free energy differences from Monte Carlo data. _Journal of Computational Physics_ 1976, _22_ , 245–268
  * Shirts and Chodera 2008 Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. _The Journal of Chemical Physics_ 2008, _129_
  * Clark et al. 2017 Clark, A. J.; Gindin, T.; Zhang, B.; Wang, L.; Abel, R.; Murret, C. S.; Xu, F.; Bao, A.; Lu, N. J.; Zhou, T. et al. Free Energy Perturbation Calculation of Relative Binding Free Energy between Broadly Neutralizing Antibodies and the gp120 Glycoprotein of HIV-1. _Journal of Molecular Biology_ 2017, _429_ , 930–947
  * Clark et al. 2019 Clark, A. J.; Negron, C.; Hauser, K.; Sun, M.; Wang, L.; Abel, R.; Friesner, R. A. Relative Binding Affinity Prediction of Charge-Changing Sequence Mutations with FEP in Protein–Protein Interfaces. _Journal of Molecular Biology_ 2019, _431_ , 1481–1493
  * Patel et al. 2021 Patel, D.; Patel, J. S.; Ytreberg, F. M. Implementing and Assessing an Alchemical Method for Calculating Protein–Protein Binding Free Energy. _Journal of Chemical Theory and Computation_ 2021, _17_ , 2457–2464, Publisher: American Chemical Society
  * La Serra et al. 2022 La Serra, M. A.; Vidossich, P.; Acquistapace, I.; Ganesan, A. K.; De Vivo, M. Alchemical Free Energy Calculations to Investigate Protein–Protein Interactions: the Case of the CDC42/PAK1 Complex. _Journal of Chemical Information and Modeling_ 2022, _62_ , 3023–3033, Publisher: American Chemical Society
  * Nandigrami et al. 2022 Nandigrami, P.; Szczepaniak, F.; Boughter, C. T.; Dehez, F.; Chipot, C.; Roux, B. Computational Assessment of Protein–Protein Binding Specificity within a Family of Synaptic Surface Receptors. _The Journal of Physical Chemistry B_ 2022, Publisher: American Chemical Society
  * Sun et al. 2021 Sun, Q.; Ramaswamy, V. S. K.; Levy, R.; Deng, N. Computational design of small molecular modulators of protein–protein interactions with a novel thermodynamic cycle: Allosteric inhibitors of HIV-1 integrase. _Protein Science_ 2021, _30_ , 438–447, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/pro.4004
  * Testa et al. 2020 Testa, A.; Hughes, S. J.; Lucas, X.; Wright, J. E.; Ciulli, A. Structure-Based Design of a Macrocyclic PROTAC. _Angewandte Chemie International Edition_ 2020, _59_ , 1727–1734, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.201914396
  * Hendrick et al. 2022 Hendrick, C. E.; Jorgensen, J. R.; Chaudhry, C.; Strambeanu, I. I.; Brazeau, J.-F.; Schiffer, J.; Shi, Z.; Venable, J. D.; Wolkenberg, S. E. Direct-to-Biology Accelerates PROTAC Synthesis and the Evaluation of Linker Effects on Permeability and Degradation. _ACS Medicinal Chemistry Letters_ 2022, _13_ , 1182–1190, Publisher: American Chemical Society
  * J. Maple et al. 2019 J. Maple, H.; Clayden, N.; Baron, A.; Stacey, C.; Felix, R. Developing degraders: principles and perspectives on design and chemical space. _MedChemComm_ 2019, _10_ , 1755–1764, Publisher: Royal Society of Chemistry
  * Ermondi et al. 2020 Ermondi, G.; Vallaro, M.; Caron, G. Degraders early developability assessment: face-to-face with molecular properties. _Drug Discovery Today_ 2020, _25_ , 1585–1591
  * Gopalsamy 2022 Gopalsamy, A. Selectivity through Targeted Protein Degradation (TPD). _Journal of Medicinal Chemistry_ 2022, _65_ , 8113–8126, Publisher: American Chemical Society
  * Towns et al. 2014 Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lathrop, S.; Lifka, D.; Peterson, G. D. et al. XSEDE: Accelerating Scientific Discovery. _Computing in Science & Engineering_ 2014, _16_ , 62–74, Conference Name: Computing in Science & Engineering