Shu-Yu Chen,Martin Zacharias
What Makes a Good Protein–Protein Interaction Stabilizer: Analysis and Application of the Dual-Binding Mechanism
Protein–protein interactions (PPIs) are essential for biological processes including immune reactions and diseases. Inhibition of PPIs by drug-like compounds is a common basis for therapeutic approaches. In many cases the flat interface of PP complexes prevents discovery of specific compound binding to cavities on one partner and PPI inhibition. However, frequently new pockets are formed at the PP interface that allow accommodation of stabilizers which is often as desirable as inhibition but a much less explored alternative strategy. Herein, we employ molecular dynamics simulations and pocket detection to investigate 18 known stabilizers and associated PP complexes. For most cases, we find that a dual-binding mechanism, a similar stabilizer interaction strength with each protein partner, is an important prerequisite for effective stabilization. A few stabilizers follow an allosteric mechanism by stabilizing the protein bound structure and/or increase the PPI indirectly. On 226 protein–protein complexes, we find in >75% of the cases interface cavities suitable for binding of drug-like compounds. We propose a computational compound identification workflow that exploits new PP interface cavities and optimizes the dual-binding mechanism and apply it to 5 PP complexes. Our study demonstrates a great potential for in silico PPI stabilizers discovery with a wide range of therapeutic applications.
What Makes a Good Protein–Protein Interaction Stabilizer: Analysis and Application of the Dual-Binding Mechanism

Protein–protein interactions (PPIs) are essential for biological processes including immune reactions and diseases. Inhibition of PPIs by drug-like compounds is a common basis for therapeutic approaches. In many cases the flat interface of PP complexes prevents discovery of specific compound binding to cavities on one partner and PPI inhibition. However, frequently new pockets are formed at the PP interface that allow accommodation of stabilizers which is often as desirable as inhibition but a much less explored alternative strategy. Herein, we employ molecular dynamics simulations and pocket detection to investigate 18 known stabilizers and associated PP complexes. For most cases, we find that a dual-binding mechanism, a similar stabilizer interaction strength with each protein partner, is an important prerequisite for effective stabilization. A few stabilizers follow an allosteric mechanism by stabilizing the protein bound structure and/or increase the PPI indirectly. On 226 protein–protein complexes, we find in >75% of the cases interface cavities suitable for binding of drug-like compounds. We propose a computational compound identification workflow that exploits new PP interface cavities and optimizes the dual-binding mechanism and apply it to 5 PP complexes. Our study demonstrates a great potential for in silico PPI stabilizers discovery with a wide range of therapeutic applications.

Complex formation of proteins plays a fundamental role in the majority of biological processes such as cell fate, immune reactions, and signal transduction and is involved in many human diseases, offering substantial therapeutic potentials for modern drug design. Approaches to protein–protein interaction (PPI) modulation are typically inhibition but more recently also stabilization of interactions. PPIs can be inhibited by allosteric effects (involving binding of compounds not necessarily at the interface) or by targeting the protein–protein interface directly such that protein partners can not form a complex. However, the design of PPI inhibitors suffers from the general feature of a typically flat and large buried interface area (>1500 Å12,13) and also often hydrophobic surface, which is frequently considered as “undruggable” for small molecules. Larger peptide-based inhibitors can be adopted as a promising strategy with low toxicity, and indeed several peptides have been approved as clinical drugs. In particular, cyclic peptides with predefined conformation and good cell permeability have received increasing attention in PPI inhibitor design, and using in silico “hot-spots” matching with cyclic peptides has been proposed for rapid rational design. For many diseases, promotion or stabilization can be more helpful than the inhibition of PPIs. Such stabilization can be achieved by small organic molecules through an allosteric (not at the interface) or a direct interface binding mechanism. The former mechanism describes the allosteric regulation of one protein partner upon stabilizer binding, which increases the binding affinity to another binding partner, while the latter mechanism describes the direct binding of the stabilizer at the PPI interface with a potentially druggable cavity for small molecules. With several advantages compared to PPI inhibitors, stabilization of PPIs has become an emerging field in recent years. Zarzycka et al. and Andrei et al. surveyed the available protein–protein-stabilizer complexes and discussed different groups of stabilizers according to the stabilized PPI type and the origin of the stabilizers, respectively. So far, only a few drug design efforts have been directed toward stabilizer design. For example, fragment-based drug discovery was used for the stabilizer design of the adapter protein 14-3-3 with p65-derived protein, the tumor suppressor protein p53, oncogenic transcription factors TAZ, and estrogen receptor-derived peptide. Sijbesma et al. performed structure-based virtual screening over nearly 6 million compounds from the Molport database and tested 13 compounds in vitro, among which 2 are validated as potent stabilizers. Tang et al. combined molecular dynamics (MD) and a molecular docking technique to investigate 22 14-3-3 protein-peptide-stabilizer complexes and concluded that a simulation time of 20–50 ns combined with molecular mechanics Poisson–Boltzmann surface area (MM/PBSA) and generalized-Born surface area (MM/GBSA) methods are useful to identify stabilizers among decoys. With the rapid advance in protein structure prediction, protein–protein docking, and protein complex structure prediction, the prediction of potential PPI stabilizers from merely the knowledge of just protein sequences can be a promising route at relatively low cost. To achieve such a goal, it is of fundamental importance to understand the functional mechanisms of the existing stabilizers. In the first part of this study, we discuss the theory of PPI stabilization and emphasize the importance of a dual-binding mechanism. Next, we analyze existing protein–protein (PP) complexes with bound stabilizers using MD simulations combined with the calculation of interaction free energies between protein partners and the stabilizer. We use the term stabilizer-induced PPI for complexes that do not form without a stabilizer (Figure 1A), whereas stabilizer-enhanced PPI includes those that form complexes already without stabilizers (Figure 1B). We find that the more potent stabilizers in most cases distribute the calculated interaction evenly between both protein partners, regardless of the total interaction free energy between the PP complex and the stabilizer. This follows the expectation based on a mathematical model for the equilibrium binding state, whereas less potent stabilizers tend to bind more strongly to one protein than the other. However, a few compounds mediate the stabilization through an allosteric mechanism by indirectly increasing the PP affinity (4 out of 18 cases). In addition, indirect positive and negative (destabilization) effects on the direct protein partner interaction were observed. Importantly, we show that 80% of the stabilizer-binding pockets can be detected in silico by direct computational probing of the PP complex (without stabilizer) crystal structure. Those pockets hidden in the PP complex structure can be revealed by running short MD simulations and subsequent pocket detection. We also demonstrate that the dual-binding mechanism can be useful to identify potential PPI stabilizers. Based on the model, we propose a protocol for PPI stabilizer discovery combining pocket probing, molecular docking, and MD simulation. We also check a large set of known PP complexes for druggable pockets and find, for 75% of the cases, interface cavities useful for stabilizer binding and apply the suggested protocol to a subset of 5 PP complexes. It indicates promising results that could be implemented in the structure-based drug design of PPI stabilizers.

Arguably, the aim of a potent PPI stabilizer is to increase the level of formation of the RLS ternary complex. In the case where R and L do not interact, e.g., stabilizer-induced PPI, the system is similar to the bivalent scaffold discussed by Yang and Hlavacek (Figure 1A). With the total concentration of receptor [R0], ligand [L0], and stabilizer [S0], the formation of the RLS complex [RLS] at the equilibrium state can be formulated aswhereϕ is the cooperative factor, [Sf] indicates the concentration of the unbound stabilizer, and KRS and KLS are the dissociation equilibrium constants of the RS and LS complex, respectively. Equation 1 describes that the RLS complex cannot be efficiently formed when [S0] is too high because of the saturation of the individual RS and LS complexes, and there exists an optimal concentration of the total stabilizer, [S0opt] = , that maximizes the RLS formation to [RLSopt] (Figure 1C). With the relation KD = eΔΔG/(RT), where KD is the dissociation equilibrium constant and ΔΔG is the Gibbs binding free energy, kB is Boltzmann constant, and T is the temperature, we can study how [RLSopt] changes with ΔΔGRS and ΔΔGLS according to eq 1. In the case where ΔΔGRS < ΔΔGLS, e.g., the stabilizer binds more strongly to the receptor, enhancing both ΔΔGRS and ΔΔGLS by a small amount improves the effective dissociation equilibrium constant KRL,eff = [R0][L0]/[RLSopt] much more efficiently than enhancing ΔΔGRS alone. This indicates that the protein that binds weaker to the stabilizer, namely, max {ΔΔGRS, ΔΔGLS, plays a more decisive role in determining the stabilization efficiency (Figure 1D). Hence, a dual-binding activity is effective for a stabilizer to efficiently induce the RLS complex formation. At the extreme case where ΔΔGRS is strong enough that all receptors are in the stabilizer-bound form, it is apparent that the formation of RLS depends solely on ligand-stabilizer interaction. Note, the dual-binding mechanism is one mechanism to stabilize PPI, and other allosteric mechanisms are discussed in the next paragraphs.

In the present study, we investigate 18 stabilizers binding to 15 different PPIs using MD simulations and pocket detection techniques. An overview of our data set is given in Table 1 and visualized in Figure 2. A stabilized PPI contains a receptor protein (R), a ligand protein (L, typically the smaller protein partner), and a stabilizer (S) and is abbreviated as an RLS complex. Similarly, a PPI without a stabilizer is abbreviated as an RL complex. The 15 PPIs with available three-dimensional (3D) structures are split into 3 data sets. Set A contains 3 stabilizer-induced PPIs where the stabilization effect of two stabilizers on the same R-L pair are experimentally measured and cocrystallized. These were selected to investigate how the experimentally measured stabilization efficiency is correlated with the structural, dynamical, and thermodynamic properties of the RLS complexes. Set B is formed by an additional 7 stabilizer-induced PPIs (structures known in complex with one stabilizer) and 5 stabilizer-enhanced PPIs from set C (for which 3D structures also in the absence of stabilizer are known). In example A1, the complex between protein 14-3-3 and the C-terminal domain of PMA2 can be induced by the dipeptide Epibestatin and (A1-a, KRL,eff = 1.8 μM) and by Pyrrolidone 1 (A1-b, KRL,eff = 80 μM). In example A2, the dimerization of the weakly dimerizing protein immunologic regulators programmed death 1 (PD-1) ligand (PD-L1, melting temperature Tm = 35.4 °C) is enhanced by small molecules BMS-202 (A2-a, Tm = 48.4 °C) and BMS-8 (A2-b, Tm = 44.8 °C). Similarly, the dimerization of bromodomain-containing protein 4 (BRD4) can be also induced by small triazolopyridazine-containing molecules compound 6 (A3-a, pKRL,eff = 8.1) and compound 6 (A3-b, pKRL,eff = 7.2) in example A3. Notably, although different binding poses of the stabilizers in A1 and A2 are observed, the RL complexes adopt similar conformations upon stabilizer binding, as indicated by the high TM-scores in Table 1. Despite differences in the interaction surface of the different stabilizers it is useful to first compare the sizes of the interfaces the stabilizers share with both protein partners. Hence, we first calculated the buried surface area (BSA) of RS and LS of the studied 18 RLS complexes (for a quick estimation of contacts between stabilizers and proteins, Figure 3A). For most stabilizers, a similar BSA is shared with both protein partners, indicating a similar number of contacts between the stabilizer and both proteins (Table 2). However, surprisingly, although A2-a and A3-a are more potent stabilizers than A2-b and A3-b, respectively, they form less BSA with R and L. Furthermore, we observed four stabilizers in A1-b, B2, B4, and B5 having only a little or no contact with the ligand, which violates the requirement for the dual-binding mechanism discussed above. To estimate the interaction more accurately, we perform MD simulations and calculated the interaction free energies ΔΔGRS, ΔΔGLS, and the total interaction free energy between the protein partners and the stabilizer ΔΔG(RL)S (Table 2 and Figure 3A) using the Molecular Mechanics Generalized Born surface area (MMGBSA) end point method (see Methods). Indeed, a comparison between the more potent and less potent stabilizers in set A shows that the weaker stabilizer-binding partner, namely, max{ΔΔGRS, ΔΔGLS}, has a more favorable calculated interaction to the more potent stabilizers (highlighted in italic font in Table 2) than to the less potent stabilizers. However, stabilizers in complexes B2, B4, and B5 do not follow this rule and also bind to the ligand-protein partner with calculated interaction free energies ΔΔGLS less than −5 kcal/mol (including also A1-b), which is the minimum binding free energy required to induce a PPI complex with KRL,eff at a μM level at [R0] = [L0] = 1 μM, according to the theoretical model (Figure 1D), implicating another possible PPI-inducing mechanism (Table 2, Figure 3B). In order to understand how these four stabilizers induce the RLS complexation without bivalent binding to both protein partners, we removed the ligand proteins and performed 200 ns MD simulations with (R + S) and without stabilizers (R-only). We found that the interface residues of the receptors, e.g., the residues in contact with the ligand protein in the RLS complexes, remain stable only in the presence of stabilizers but become much more mobile when the stabilizers are removed, indicated by increased interface root-mean-square deviation (iRMSD, Figure 3C–D). Hence, in these cases, the stabilizers bind to one protein partner (receptor) and stabilize the bound conformation, and in turn, allosterically stabilize PP complex formation. We next calculated the direct interaction between receptor and ligand with the stabilizer and without (ΔΔGRL) the stabilizer (see Methods). Figure 3E shows that the calculated RL interaction is weakened in most stabilizer-induced PPIs when the stabilizers are removed but is stronger or similar in most stabilizer-enhanced PPIs without stabilizers. In order to understand this effect, we performed residue-wise energy decomposition analysis on the RL complexes simulations and identified the favorable RL contacts (indicated by negative binding energies) and unfavorable RL contacts (indicated by positive binding energies, SI Appendix, Figures S1 and S2). Interestingly, the unfavorable RL contacts in the stabilizer-induced PPI complexes are often formed by the residues responsible for stabilizer binding (SI Appendix, Figures S1 and S2). Thus, in addition to the dual-binding mechanism, stabilizers can also contribute to more favorable binding by shielding the unfavorable RL interactions. This may explain why the conformations of receptor and ligand in set A tend to form a similar conformation upon stabilizer binding where favorable RL contacts and druggable pockets near the unfavorable RL contacts can be formed simultaneously (Table 1). When stabilizers are removed from these RLS complexes, these unfavorable contacts are exposed to their binding partner, causing interface distortions. Such an effect was observed in complexes A3, B6, and B7 within 50 ns simulations (Figure 3F, SI Appendix, Figure S3, and SI Appendix, Figure S4). In contrast, the interfaces of the stabilizer-enhanced PPIs remained stable without stabilizers (Figure S5). The binding of these stabilizers involves many favorable RL contacts (Figure 3G and SI Appendix, Figure S6) and is therefore more likely to disturb the direct RL interactions.

In order to discover potential PPI stabilizers in the scheme of structure-based drug design (SBDD), it is crucial to identify useful potential ligand-binding pockets. Structural comparison has shown that PPI stabilizers can bind into an interface pocket of an RL complex without deforming the conformation (Table 1), offering an opportunity for virtual screening. We used Fpocket4.0 to probe druggable pockets of the natural RL complexes and detected at least 10 pockets in each RL complex (Figure 4A). Although the dual-binding mechanism already excludes the noninterface pockets, performing virtual screening through all possible interface pockets is still computationally expensive, and therefore, a reliable approach to rank the interface pockets is desired to narrow down the sampled phase space. To examine whether the scores provided by Fpocket, e.g., the drug score and pocket score, can accurately distinguish the stabilizer-binding pocket from other interface pockets, we calculated the ligand coverage fraction of each pocket (see Methods). Besides the C4 complex, all stabilizer-binding pockets are fully or partially detectable from the RL complex and well ranked by both pocket score and drug score, indicating that it is a reliable approach to identify potential stabilizer-binding pockets solely on the basis of the RL structure (Figure 4B). Further the use of MDpocket to track the dynamics of stabilizer-binding pockets from the simulation (see Methods) reveals that all pockets exhibit a breathing motion with a considerable range of volume sampled in the simulations (Figure 4C and time series shown in SI Appendix, Figure S7). Importantly, the stabilizer-binding pocket hidden in the RL crystal structure of the C4 complex occasionally expands to a volume large enough to accommodate the stabilizer with a rapid fluctuation in pocket volume (Figure 4C and SI Appendix, Figure S7). Upon removal of the stabilizer from the RLS complexes, the volume of the stabilizer-binding pocket showed in a few cases little changes during simulations. However, in many cases, a reduced pocket size distribution sampling frequently pocket volumes smaller than the ligand volume was observed (Figure 4D). In the A3, B6, and B7 cases, the removal of the stabilizer resulted in a distortion of the RL interface such that the stabilizer binding pocket became completely undetectable (SI Appendix, Figures S3, S4, S8, and S9). Note, interface distortion and RL complex dissociation are expected in all stabilizer-induced PPI if simulations were extended to much longer (e.g., microsecond) time scales due to the very weak interaction in these cases.

Our systematic study of 18 PPI stabilizers reveals PPI stabilization mechanisms and the dynamics of the stabilizer-binding pockets. In particular, the dual-binding mechanism derived from the starting theoretical model and supported by experimental results shows that the binding affinity between the stabilizer and the weakly interacting protein partner plays a decisive role in the RLS complexation process. Hence, it may provide a better ranking than the conventional structure-based drug design (SBDD) workflow, which usually only considers the total calculated interaction energy between the compound and the target of interest. We tested the dual-binding ranking approach on the 14-3-3/ChREBP complex with 13 compounds preselected by Sijbesma et al. through a virtual screening approach (SI Appendix, Figure S10A). We performed molecular docking, MD simulations, and MM/GBSA calculations to check the dual-binding mechanism of these compounds (see Methods, Supplementary Figure 10B). For all compounds, a stronger calculated interaction with the receptor compared with the ligand was observed (Figure 5A). According to the proposed dual-binding mechanism, the interaction of the compound and the weaker ligand partner then plays a more decisive role in stabilization efficiency. Indeed, an improvement in ranking performance (compared to the experiment) was found when the compounds are ranked using the binding affinity to the ligand instead of the receptor or vs total binding energy (Figure 5A). To this end, we propose an in silico protocol adapted for PPI stabilizer discovery (Figure 5B). A 3D structure of the desired RL pair needs to be available, e.g., experimental structure or based on protein–protein complex modeling. Next, Fpocket can be used to probe and rank druggable pockets at the PPI interface. One can further assess the stability of the pocket or explore hidden pockets by performing MD simulations. Once a suitable binding pocket is identified, virtual screening can be conducted with molecular docking programs to screen through a large chemical database (e.g., following the protocol described in the Methods). Compounds with high docking scores can then be submitted to MD simulations to study the binding stability and interaction free energy. At this step, it is important to check if the compound-protein interaction is of a similar magnitude with both protein partners (desired to optimize the stabilizing effect according to the dual binding mechanism). Finally, compounds that fit the criteria can be submitted to experimental validation and lead optimization. As a demonstration, we implemented the protocol to a subset of protein docking benchmark set 5.0, that contains 226 high-quality protein complexes. Pocket detection revealed that among 226 tested complexes, 223 complexes (98.7%) have at least one pocket consisting of 5 atoms from each side of the interface, 173 complexes (76.5%) have at least 10 atoms from each side of the interface, and 57 complexes (25.2%) have at least 20 atoms from each side of the interface (SI Appendix, Figure S11A–B). We further extracted the best-scored interface pockets composed of at least 5 atoms from each protein, and several pockets were identified, some with pocket scores and drug scores higher than the known stabilizer-binding pockets discussed in the previous sections (Figure 4B, SI Appendix, Figure S11C–D). Among these pockets, we chose 5 interface pockets with high drug scores, including PDBID 1AK4, 1Y64, 2GAF, 3H11, and 4FZA, for virtual screening (Figure 5C). We performed molecular docking on these five interface pockets with 1 million compounds selected from the ZINC20 database (see Methods). Besides the pocket predicted with the poorest drug score (score = 0.812, PDBID 3H11), all pockets can accommodate over 75% of the screening compounds with a docking score better than −6 kcal/mol, indicating good druggability preselected with drug scores (Figure 5D). The top-ranked 10 docked compounds of each protein complex and their docking score are listed in SI Appendix, Figures S12–S16 and were submitted to MD simulations for interaction free energy calculation (SI Appendix, Tables S1–S5). As demonstrated in the present study, an ideal PPI stabilizer usually possesses a strong binding affinity to the weaker stabilizer-binding partner, indicated by a low max{ΔΔGRS, ΔΔGLS}, and a similar binding affinity to both proteins, indicated by a low |ΔΔGRS – ΔΔGLS|. Among the 50 selected compounds, 34 compounds (68%) possess a max{ΔΔGRS, ΔΔGLS} lower than −15 kcal/mol, and a |ΔΔGRS – ΔΔGLS| smaller than 10 kcal/mol (Figure 5E). In particular, the best three compounds at the PP interface of 4FZA (compound 5, compound 8, and compound 9, SI Appendix, Table S5) all exhibited a max{ΔΔGRS, ΔΔGLS} lower than −24 kcal/mol, and |ΔΔGRS – ΔΔGLS| smaller than 4 kcal/mol, and are expected to be promising stabilizers. Together with the dual-binding mechanism of stabilizers, our results provide not only a prediction of potential PPI stabilizers but also a guide for further lead optimization, for example, by improving the RS interaction of compound 3 in PDBID 1AK4 or by improving the LS interaction of compound 4 in PDBID 3H11 (SI Appendix, Table S1 and Table S4).

To facilitate the development of the PPI stabilizer discovery, it is essential to understand how these compounds stabilize PP complexes. In the present study, we have case-studied the binding behavior of 18 stabilizers on their corresponding PP complexes using MD simulations, interaction free energy calculations, and pocket detection techniques. Stabilizers can facilitate PPI complex formation through the dual-binding mechanism and/or conformational stabilization of the bound receptor (or ligand) structure. The more common optimal dual-binding mechanism requires the stabilizer to interact approximately equally with the protein partners. The effective stabilization depends predominantly on the partner that binds more weakly to the stabilizer. Furthermore, stabilizers can shield the unfavorable contacts between a noninteracting RL pair while not affecting the favorable contacts. We also found that a fraction of stabilizers that do not fulfill the dual-binding mechanism may also support PPI complexation by stabilizing the interface of one protein to a conformation ready for partner protein binding. The latter more complex allosteric mechanism could, in principle, be detected by extensive MD simulations of a protein partner in the presence and absence of the stabilizer. However, the result of such simulations can be sensitive to force field artifacts, and due to the computational demand, it may not be useful for investigating many compounds. Hence, in practical design efforts, it might be most promising to focus on identifying putative stabilizer compounds that follow a dual-binding mechanism. Indeed, our analysis of the stabilizer-binding pockets in experimental structures without stabilizers demonstrated that most stabilizer-binding pockets can be readily detected. Those pockets hidden in the RL complexes can show up in short MD simulations. Furthermore, interface pockets that can accommodate ligands, allowing stabilization by the dual-binding mechanism, could be detected in the majority of cases in a large set of known PP complex structures. Based on our analysis, we demonstrated that the dual-binding mechanism can be useful to identify potent stabilizers and propose a protocol for PPI stabilizer discovery. We applied the protocol to 5 test cases and can suggest potential PPI stabilizers for five protein–protein complexes. Combining the pocket-scoring technique and dual-binding criteria, our SBDD workflow can efficiently identify druggable interface pockets and estimate the stabilizing efficiency in silico. Ideally, with advances in structural prediction and protein–protein docking techniques, one can design molecules that induce the complexation of selected proteins. For a given interaction geometry, our protocol could also be useful for the identification of such compounds that induce a PPI or support a very weak association. Our work is the first computational study systematically characterizing the binding behavior across different kinds of PPI stabilizers. We anticipate that our study will facilitate future developments for the discovery of PPI stabilizers.

In this study, in total, 23 protein complex structures were used from the RCSB PDB database, including 18 RLS complexes and 5 RL complexes (Table 1). Due to missing side chains and short backbone segments in some proteins (and some mutated residues), the comparative modeling software MODELER was used with the wild-type sequences listed in Table 1. This includes the restoration of the mutations N755S in PDBID 3bbr, T686A in PDBID 3b6q, L22Q, T51D, V52I in the ligand of PDBID 3o98, S13A, T30D, V31I in the ligand of PDBID 3m50 and 3m51, C6S in PDBID 3u15, and L9E and 3vbg. The RL complexes of sets A and B are generated by stripping the stabilizers from the corresponding RLS complex. Because of the high structural similarity, A1-a, A2-a, and A3-a were taken to represent the RL complexes of A1, A2, and A3, respectively. Besides the 14-3-3 protein and its ligands, where the terminal residues are as well interacting with each other and the stabilizers, the N-terminus and C-terminus are capped with the neutral capping group ACE and NME with Ambertools 20, respectively. To maintain the conformations of the complex, two calcium ions are kept in the S100A4 homodimer and CaM/CaMBD2a complex, four glutamic acids are kept in the iGluR2 homodimer, and a zinc ion is kept in the CK1α/CRL4. Also, disulfide bonds in PD1L/PD1L and lambda-6A/lambda-6A complexes are preserved. Each system was solvated in a box 12.5 Å extended from the nearest atom of the solute with OPC 4-point water molecules and 0.15 M sodium chloride. The protonation state of the stabilizers is determined with Open Babel at pH = 7.4. The atomic interactions are described by the ff19SB force field for protein, OPC force field for water, the TIP4PEW force field for ions, and GAFF2 for the stabilizers. The same procedure was also implemented to prepare the MD simulation of the 5 example complexes from protein docking benchmark data set 5.0. Five replicas of 50 ns trajectories of each RL and RLS complex and 5 replicas of 200 ns long trajectories of each RS complex and receptor-only system of PDBID 3m51, 2o98, 4mdk, and 1kkq are generated with MD simulation with the CUDA-accelerated version of PMEMD from AMBER20 package with randomly assigned initial velocities. Each simulation box underwent the minimization process of maximal 35,000 steps and was equilibrated for 75 ps under constant volume and 225 ps under constant pressure. Gradually decreasing positional restraints on protein atoms are applied during the minimization and equilibrium processes. The temperature was kept at 310 K using a Langevin thermostat, and the pressure was kept at 1 bar with a Berendsen barostat. The nonbonded cutoff distance was set to 9 Å, and the simulation was integrated with a time step of 2 fs with the SHAKE algorithm. Interface residues are defined by the residues within 5 Å of the partner protein calculated using VMD-PYTHON. RMSD and iRMSD are calculated with the starting structure as the reference. RMSF, RMSD, and iRMSD are calculated using CPPTRAJ. The last 30 ns, in a total of 75 frames, of each 50 ns simulation in the RL and RLS simulation are taken to calculate the interaction free energies between receptor, ligand, and stabilizer using molecular mechanics coupled with the generalized-Born surface area (MM/GBSA) method with an internal dielectric constant of 1 and external dielectric constant of 80. The effective Born radii are calculated with the GBOBCII model (igb = 5) in ref (65). Pocket detection was performed with Fpocket4.0 on the crystal structure of the RL complexes in set C. Each pocket consists of polar and apolar probes, as shown in Figure 5C in red and gray spheres, respectively. An atom and a pocket probe are considered in contact if the distance between the atom and the pocket probe is smaller than the sum of their radii with the radii of atoms set to H: 1.2 Å, C: 1.7 Å, N: 1.55 Å, O: 1.52 Å, F: 1.47 Å, B: 1.92 Å, P: 1.8 Å, S: 1.8 Å, Cl: 0.2 Å. L. The ligand coverage fraction is calculated as , where Nligand ∩ Npocket is the number of ligand atoms contacting the pocket probes, as illustrated in Figure S17. The dynamics of the stabilizer binding pocket is analyzed using MDtraj by first creating probes within 2 Å from the known stabilizers, as shown in the SI Appendix, Figure S18. Re-evaluation of the 13 compounds performed by Sijbesma et al. was generated with Open Babel and docked with Autodock vina. The docking procedure was performed on the 3D structure from PDB 5F74. The preparation of 5*50 ns trajectories and the binding energy calculation were performed with the same procedure as mentioned above. Considering many stabilizers in the present study possess molecular weights between 375 and 425 Da, and LogP between 2 and 4.5 (SI Appendix, Figure S19), we fetched 122,847,475 compounds fulfilling such molecular properties from the ZINC20 database with the standard highest reactivity. We further eliminated compounds with complex structures, which might cause difficulties for force field parametrization, by filtering the compounds with BertzCT complexity > 700 using RDkit, leading to 15,072.167, from which we randomly selected 1 million compounds to perform molecular docking with AutoDock-GPU.