Characterization of Promiscuous Binding of Phosphor Ligands to Breast-Cancer-Gene 1 (BRCA1) C-Terminal (BRCT): Molecular Dynamics, Free Energy, Entropy and Inhibitor Design

Inhibition of the protein-protein interaction (PPI) mediated by breast-cancer-gene 1 C-terminal (BRCT) is an attractive strategy to sensitize breast and ovarian cancers to chemotherapeutic agents that induce DNA damage. Such inhibitors could also be used for studies to understand the role of this PPI in DNA damage response. However, design of BRCT inhibitors is challenging because of the inherent flexibility associated with this domain. Several studies identified short phosphopeptides as tight BRCT binders. Here we investigated the thermodynamic properties of 18 phosphopeptides or peptide with phosphate mimic and three compounds with phosphate groups binding to BRCT to understand promiscuous molecular recognition and guide inhibitor design. We performed molecular dynamics (MD) simulations to investigate the interactions between inhibitors and BRCT and their dynamic behavior in the free and bound states. MD simulations revealed the key role of loops in altering the shape and size of the binding site to fit various ligands. The mining minima (M2) method was used for calculating binding free energy to explore the driving forces and the fine balance between configuration entropy loss and enthalpy gain. We designed a rigidified ligand, which showed unfavorable experimental binding affinity due to weakened enthalpy. This was because it lacked the ability to rearrange itself upon binding. Investigation of another phosphate group containing compound, C1, suggested that the entropy loss can be reduced by preventing significant narrowing of the energy well and introducing multiple new compound conformations in the bound states. From our computations, we designed an analog of C1 that introduced new intermolecular interactions to strengthen attractions while maintaining small entropic penalty. This study shows that flexible compounds do not always encounter larger entropy penalty, compared with other more rigid binders, and highlights a new strategy for inhibitor design.


Introduction
The tandem~100-amino acid repeats of breast-cancer-gene 1 (BRCA1) C-terminal (BRCT) are known to bind to phosphorylated proteins which are important for a number of tumor suppressor functions, which include, DNA repair, cell-cycle checkpoint, and transcription regulation [1][2][3][4]. The BRCT repeats recognize and bind phosphorylated protein partners such as BRCT uses to achieve promiscuous binding and the interaction energy of the ligand-BRCT. The MD simulations illustrated the molecular flexibility in the free and bound states for BRCT (BRCA1) and ligands. We analyzed loop movements and the population of dihedral rotations of backbone and side-chains. Conformations from MD simulations were used as initial structures for thorough conformational search and free energy calculations with the M2 method, to reveal the contribution of configuration entropy and enthalpy to ligand binding affinities. We focused on how to optimize the balance between enthalpy gain and entropy loss. Using an accepted practice in ligand design, we synthesized a ligand that incorporates a benzene ring to possibly constrain its conformation. Upon ligand binding, changes of each energy term, conformations, rotameric state, and configurational entropy were evaluated by both MD and M2 tools; and the findings were used to suggest new inhibitors. Table 1 lists 14 short peptides (P1-P14, among which P11 contains a phosphate mimic and others contain phosphorylated amino acids) [24], one compound (C1) [27], one new compound (N1), one designed compound (D1) and 4 long phosphopeptides (L1-L4) that bind to the BRCT domain: pS and pT is phosphorylated amino acid serine and threonine, respectively, and γcE is γ-carboxyglutamate, which is chosen to mimic pS interaction as a non-phosphorylated peptide binder.

MD simulations
We ran MD simulations on BRCT-ligand complexes, free ligands and free protein, and the PDB IDs used as initial structures to perform MD simulations were listed in S1 Table. The initial bound conformation of all tetrapeptides was generated by superimposing the backbone atoms of -pSXXF-within phosphorylated BACH1 peptide ISRSTpSPTFNKQ in the C-terminal domain of the BRCA1 protein (PDB code 1T29) [46]. Besides 1T29, we included the other three BRCT domain structures in complex with long phosphopeptides from CtIP, ACC1 proteins and library screening, with PDB IDs 1Y98 (PTRVSpSPVFGA), 3COJ (PQpSPTFPEAG) and 1T2V (AAYDIpSQVFPFA), respectively, for promiscuous molecular recognition study [47][48][49]. The initial structure of the bound conformation of C1, N1 and D1, where no available crystal structures, were from docking with Autodock tools 1.5.6 [50,51] and then further checked manually by ensuring important interactions hold. Notably, Autodock was used for only the three ligands that did not have co-crystal structures with BRCT. The docking method used the Lamarckian genetic algorithm, which fixed the protein and allowed the ligand to move around in the docking box. The partial charges of ligands were calculated by using the Vcharge program [52]. The Autodock scoring function is a subset of the AMBER force field that treats molecules using the united atom model. Autogrid version 4.0 was used to create affinity grids with 0.375 Å spacing in 19.5 x 11.25 x 11.25 Å 3 space at binding site. The final docking result was obtained by 10 runs of simulation with 2.5 million rounds of energy evaluation in each run. Ligand conformations with the lowest docked energies and reasonable conformation (pSer forms hydrogen bonds with S1655, G1656 and K1702, and the P+3 Phe locates in the hydrophobic packet formed by M1775, N1774 and F1704) were further analyzed. We selected two initial conformations with similar low energy computed by Autodock for ligands C1 and D1, and N1 has one initial conformation (S1 Fig).
We performed MD simulations on an apo BRCT domain, 21 complexes, and 21 free ligands to study the dynamic nature of a given system. The standard simulation package, Amber14 [53] with the Amber 99SB force field [54][55][56][57], was used. For pSer and pThr, we used the force field reported by Homeyer et al [58]. Amber atom types were manually assigned to non-standard amino acid and functional groups of the ligands C1, N1 and D1. Each system was set up as follows. First, we minimized the hydrogen, side-chain and whole system for 500, 5 000 and 5 000 steps, respectively; then the systems were solvated in a rectangular box of a 12-Å explicit TIP3P water model by the tleap program in Amber14. Each system contains about 50 000 atoms. Counter ions Na + were added to keep the whole system neutral, and particle mesh Ewald was used to consider long-range electrostatic interactions [59]. Before equilibration, we ran energy minimization of 10 000 and 20 000 steps for the waters and system, respectively; next, we ran equilibrium of solvent molecules for 40 ps. Then the systems were gradually heated from 250 K for 20 ps, 275 K for 20 ps, to 300 K for 160 ps. We saved a frame every 1 ps with a time step of 2 fs in the isothermic−isobaric (NPT) ensemble. The Langevin thermostat with a damping constant of 2 ps −1 was used to maintain a temperature of 300 K, and the The major binding residues, pSer and Phe (P+3), are in bold. The relative binding free energy for ligand X (X = P2-P14, C1, N1) to ligand P1 is approximated using the half maximal inhibitory concentration IC50 as ΔΔG exp = RT ln IC50(X)/IC50(P1) based on equation ΔG = RT ln K d = RT ln(IC50 + 0.5C enzyme ) % RT ln IC50 [44,45]. Binding free energies for L1-L4 are calculated through equation ΔG exp = RT ln (K d ). a IC50 values of P1-P14 were taken from ref [24]. IC50 values of C1 was taken from ref [27]. b K d values of L1-L4 were taken from ref [46], ref [47], ref [48] and ref [49], respectively. hybrid Noseé−Hoover Langevin piston method was used to control the pressure at 1 atm. We also used the SHAKE procedure to constrain hydrogen atoms during MD simulations [60]. Finally, all production runs were performed for 100 ns at 300 K. To ensure that all simulations reached stable energy fluctuations, we considered only trajectories during 20−100 ns for postanalysis.

M2 method
The second-generation mining minima method, M2, calculates the standard free energy of binding by computing the free energy of the free BRCT (G°B RCT ), ligand (G°l igand ), and ligand-BRCT complex (G°c omp ).
M2 uses the classical formulation of the partition function for calculating free energy G°.
where U is potential energy, W is the solvation free energy and Z i is the local configuration integral from distinct energy wells. The external degrees of freedom were integrated out and C°p rovides a correction to the standard state, and r_int indicates the variables of the internal bond-angle-torsion coordinates. Formally, the configuration integral must be determined over all spaces along the remaining internal degrees of freedom. M2 approximates this configuration integral by using the concept of considering local energy minima only [61,62]. Therefore, the M2 approach replaces the configurational integral over all spaces with a sum over separate local configurational integrals (Z i ) associated with the low energy minima of the system. Determining Z i allows for the probability to be associated with each energy well, which in turn, allows for determining a Boltzmann averaged energy <U+W>, which is then subtracted from the total free energy to give the system configurational entropy, useful when analyzing and interpreting predicted binding affinities.
Note that the configurational entropy S°c onfig includes both a conformational part, which reflects the number of energy wells (conformations), and a vibrational part, which reflects the average width of the energy wells. The solvent entropy is included in the solvation free energy, W. Therefore, the computed configurational entropy changes cannot be directly compared with experimentally measured entropy changes, which contain both configurational and solvent entropy.
In brief, M2 contains two parts: 1) an aggressive conformational search for distinct lowenergy wells, with repeats detected and removed; and 2) an enhanced harmonic approximation for computing the configuration integral Z i of each well i. Each distinct conformation is energy minimized, first by conjugate gradient method and then by Newton-Raphson method. Both parts involve the Hessian matrix with respect to bond-angle-torsion coordinates, and our harmonic approximation accounts for anharmonicity of eigenvectors of the Hessian matrix with eigenvalues < 2 kcal/mol/Å or 2 kcal/mol/rad. The correlation between different degrees of freedom (e.g., multiple dihedrals may rotate in concert or move with ligand translation/rotation) is captured in the Hessian matrix. We used the VM2 package for the calculation [63][64][65] and performed three iterations for each ligand and 3 to 10 iterations for the free BRCT and the complexes until the cumulated free energy converged (S2 Fig). To reduce the computational cost, only parts of BRCT were flexible, called the "live set" (Fig 3), which are residues within 7 Å of a long peptide ISRSTpSPTFNKQ in complex with BRCT (PDB code 1T29). The rigid set, called the "real set", contained the residues within 5 Å of the live set. Other atoms not included in these two sets were not considered during the M2 calculations. All ligands were completely flexible and can freely translate and rotate within the binding site, and the same rigid and flexible parts of BRCT were applied to all systems.

Post-MD analysis: Identifying rotamer states and MM/PBSA calculations
To compare the conformational changes of a molecular system between its free and bound states, we analyzed the selected ligand and BRCT dihedral angles during MD simulations and M2 calculations. Dihedral angles were measured by using T-analyst [66], which can detect the angle population to find discontinuity in a dihedral distribution such as one energy well splitting into two wells near -180°and +180°. A shifted angle by adding or subtracting 360°is then applied to illustrate proper rotamer states. The population of each dihedral was then plotted by using Matlab with a histogram of 144 bins ranging from -360°to +360°to ensure coverage of all rotamer states after angle shifting. When analyzing the rotameric states, because the analysis does not need more than 1000 data points [66], we used trajectories with a smaller file size that a frame was saved every 100 ps (1000 frames) for each 100 ns MD run.
We used the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA)-type postprocessing method to compute ligand-BRCT intermolecular interactions during MD simulations [67][68][69][70][71][72][73][74][75][76]. The interaction energy, Δ(U+W) associated with BRCT and a ligand is computed by Δ(U+W) = <E complex >-<E bound BRCT >-<E bound ligand >. The bracket <E> denotes the average energy computed from a given MD trajectory and the energy terms include a valence term (bond, angle and dihedral), van der Waals (U VDW ), Coulombic (U Coul ), solvation free energy computed by the Possion-Boltzmann equation (W PB ) and by cavity/surface area (W NP ). The dielectric constants of the interior and exterior protein were set to 1 and 80, respectively. The valence term was canceled because of the single trajectory approach.

Ligand N1 synthesis, purification, and determination of IC50
Peptide synthesis. The peptide with the modified amino acid was synthesized by using standard Fmoc chemistry following previously reported methods [24,77,78]. The peptide was purified by preparative LC to >95% as assessed by HPLC and characterized by mass spectrometry (S3 Fig).
Protein expression and purification. The plasmid construct (pAM15, gift from Luc Gadraeu, University De Sherbrook) encoding six his-tagged BRCT domains of BRCA-1 (amino acids 1646-1859) was used to transform BL21(DE3) RIL (Stratagene). Protein expression was induced by 1mM IPTG and the recombinant protein was purified by nickel affinity chromatography (Qiagen). Homogeneity of the purified protein preparation was assessed by SDS-PAGE and concentration estimated by BCA method (Pierce).
Fluorescence polarization assay. The peptide was evaluated in a BRCT assay following previously reported methods [79][80][81], representative dose-response curves from our previous Fluorescence polarization assay study was shown in S4 Fig. It was carried out in a 384-well low volume corning plate. The polarization and fluorescence were measured on a Spectramax M5 (molecular devices) plate reader. The peptide was titrated into a mixture of BRCT(BRCA1) (1000 nM) and Fluorescently labeled peptide Flu-βA-pSPTF-CONH 2 (100 nM) where βA is beta-alanine. The IC50 value was calculated by using SigmaPlot. Unfortunately, N1 peptide was inactive even at 1000 μM (1 mM) concentration and all we got was a flat line.

Results and Discussion
We first applied MD simulations and post-MD analysis for the peptides (P1-P14, L1-L4) and compound C1 to study the fluctuations in various complexes, followed by more rigorous free energy calculations with the M2 method for short peptides (P1-P14) and compound C1 to illustrate detailed energetic and entropic changes upon ligand binding. The new ligand N1 based on consensus ideas that impose structure constraints, was examined experimentally and computationally. Based on our results, compound D1 was derived from the tight binder C1.

Conformational flexibility of the molecular systems
One unique feature of promiscuous protein systems such as BRCT is to bind to various ligands with significantly different size and shape by using the same binding interface. BRCT needs to provide adequate conformational isomers to recognize these ligands, which involves both sidechain rotation and additional plasticity provided by the backbone. As what shown in crystal structures of BRCT, the relatively rigid alpha helix and beta sheets hold the overall geometry. The variety of side-chains of residues in loops (β3-α2 connection loop, β1'-α1' connection loop and linker between N-terminal and C-terminal) creates a binding surface for ligand recognition except for the reserved binding region for the phosphate group [46]. The backbone nitrogen of G1656 and side-chain of S1655 of the β1 sheet and K1702 of the α2 helix form at least three stable hydrogen bonds with the phosphate group and also orient a ligand in the binding site (Fig 1). Notably, the pocket reserved to bind the phosphate group is located between a structurally rigid region constructed by a helix and a sheet. In contrast, the hydrophobic pocket for the P+3 phenylalanine is built by M1775 and N1774 of the β1'-α1' connection loop and F1704 of the α2 helix, with the β1'-α1' connection loop providing a certain flexibility for peptide binding [48].
To study the flexible regions in the binding pocket of BRCT, we measured the root mean square fluctuation (RMSF) of C α and the standard deviation of phi and psi angles of residues in the BRCT backbone within 7 Å of 18 peptides (P1-P14, L1-L4) and compound C1. The RMSF in S5 Fig shows that residues contacting with a ligand generally have smaller fluctuations and residues without contact with a ligand generally have larger fluctuations, Except for P13, where the middle two proline residues of tetrapeptides do not form optimized contacts with BRCT. Although RMSF plot suggested that residues contacting with a ligand have small fluctuations in the Cartesian space, the standard deviation of phi and psi angles in Fig 4(A) shows that the backbone dihedral angle can still rotate considerably. As illustrated in Fig 4(B), the most flexible region in the center part of the binding pocket, which directly contacts with the middle two residues of a pSer-X-X-Phe peptide and middle atoms of compound C1. Utilizing the flexible loop region allows for the polar residues E1698 and R1699 of the β3-α2 connection loop to form a hydrogen bond with backbone atoms of the phosphopeptides and also accommodate ligands with different shapes. For example, the standard deviations for E1698, R1699 and T1700 were especially large when BRCT bound to P13 and C1, followed by concerted motions of N1742 and G1743 in the linker region. Although P13 still can fit into the binding cavity, the two proline residues limit the arrangement of both molecules to optimize the intermolecular interactions. In contrast, C1 was flexible and adopted multiple bound conformations to strengthen its binding affinity, as discussed in the following sections. For the long peptides, F1772, T1773 from the β1'-α1' connection loop and D1692, A1693 from the β3-α2 connection loop fluctuate to adjust the size of the binding cavity. The size change of binding site agrees with our previous molecular dynamics study, where the size of cavity can be characterized by two angles E1698-A1752-E1836 and S1655-A1752-N1774, which can have difference of 10°u pon binding of different peptides (S6 Fig) [25]. In summary, BRCT uses the power of loops to alter the shape and size of the binding site to fit various ligands, combined with a rigid region designed to form stable hydrogen bonds with the phosphate group.

Ligand binding modes and intermolecular interactions computed by MM/ PBSA calculations
Because the BRCT domain has a highly adaptable binding pocket, we hypothesized that some ligands may feature diverse binding modes. We therefore examined the ligand binding modes and the rotamer of each rotatable bond for every ligand to discover their differences between the free and bound states. For all peptides P1-P14, only one major bound conformation was observed: pSer forms hydrogen bonds with S1655, G1656 and K1702 and the P+3 Phe locates in the hydrophobic packet (Fig 1). Interestingly, compound C1 can establish multiple bound conformations in the binding site by fitting either a benzene ring into the hydrophobic pocket and an indole ring into a cluster of residues G1656, L1657, T1658 of the β1-α1 connection loop and K1690 of the β3-α2 connection loop, and vice versa (Fig 5(C) and 5(B)). C1 can also bind to BRCT with its folded form, whereby two rings form a T-shape stacking interaction (Fig 5(A)).   ligands and BRCT are about the same, which agrees with experiments finding that ΔΔG exp is within 3 kcal/mol (Table 2).
To understand why or why not a ligand loses the rotamers after binding, we clustered conformations of the free peptides and ligands and compared them with those in the bound complexes. For the peptides and C1, they generally have two distinct conformations in the free state, folded and extended, which except for P13 (Ac-pSPPF-NH 2 ), can switch back and forth in MD simulations of free ligands (Fig 7). However, the bound peptides are locked to only the extended form, which results in reduced rotamers in side-chains and also backbone φ and ψ angles (Figs 6 and S7). To test the robustness of MD simulation on rotameric states analysis, we ran and analyzed another MD run with different initial conformations for several ligands. The simulated  5).

Binding free energies with M2 method
To gain insights into the mechanism of binding, we needed thorough sampling and accurate ligand binding free energy calculations that included both enthalpic and configurational entropic contributions for molecular recognition. Although MM/PBSA calculations provide valuable information for intermolecular interactions, our calculations based on 100-ns MD simulations may have missed some important conformations, and contributions from changing configurational entropy and molecular conformations are neglected in Table 2. In addition, because of different non-polar solvation models and use of a real set in M2 for energy calculations (Fig 3), the values of non-polar and polar interaction energies, ΔE NP and ΔE Polar , from MM/PBSA and M2 cannot be compared directly. We therefore computed ligand-binding free energy with the M2 method, which involved an aggressive conformational search engine to locate local energy minima and a rigorous modified harmonic approximation approach to compute free energy for each minimum found. Table 3 and Fig 8 show that the computed related binding free energy, ΔΔG calc , was in good agreement with experimental values, which validated the method as well. Because M2 uses accumulated energy which is different from dynamics-based method, it does not have fluctuated energy. Eq 2 shows that when a low energy minimum is found by M2 conformational search and added to the accumulated energy, the computed free energy G°drops. Search and The binding free energy was computed by ΔG cal = G complex -G free protein -G free ligand . Each decomposed energy is obtained by ΔE cal = E complex -E free protein -E free computation continue until the accumulated free energy is converged. Here we calculated error interval for y-intercept of linear regression line in Fig 8 [83 -85]. Part of the variance comes from experimental noise, which is typically about 0.3-0.5 kcal/mol for accurate binding free energy measurements [86]. If the binding free energies of two ligands are measured independently in experiments, then experimental relative binding free energies between the two ligands would have error around 0.4-0.7 kcal/mol. Therefore, the errors for free energy calculation method versus experimental data can only be larger than experimental noise of 0.4-0.7 kcal/ mol, indicated by the range of y-intercept of linear regression line (~3 kcal/mol), and experimental noise is expected to be a at least 13% of the total observed error. With agreement of early studies on ligand-protein binding, the strong Coulombic attraction is largely compensated by the solvation free energy, and the vdW attraction is the major driving force for ligand binding [32,64]. Moreover, peptides with large non-polar residues at the P+2 position, such as P2, P3, P5, P8, P10 and P12, generally have stronger vdW interaction (Table 3). Although M2 revealed more bound conformations for the complex from various combinations of side-chain rotations, the major binding mode of BRCT-pSXXF is the same as that obtained by MD sampling, whereby the phosphate group forms hydrogen bonds with S1655, G1656 and K1702, and P+3 Phe or Tyr locates in the hydrophobic pocket (S10 Fig). M2 also revealed more conformations for free ligands, including the folded and extended forms, and their computed conformational free energies are similar. Therefore, the folded and extended conformations may have similar population in the free ligands. It is not surprising to observe the enthalpy Δ<U+W> and configuration entropy-TΔ<S> compensation for tight binders; however, the outlier C1 is particularly of interest (Fig 9). Although C1 forms a moderate enthalpy attraction with BRCT,~-38 kcal/mol, which is similar to that for peptides P4-P9, with the remarkable~2-4 kcal/mol small configuration entropy loss, C1 outperforms other peptides (Table 3). Compared with peptides, some rotamers of BRCT and C1 can gain new rotameric states rather than losing them, and the vibrational entropy loss is smaller than that for P4-P9, as seen from the change in width of M2 histogram peaks that correspond to the width of energy wells (Figs 6 and S7). Upon ligand binding, M2 histogram peaks for P1-P14 become narrower, whereas C1 has the same or even wider peaks. In S2 Table, we list the number of complex, ligand and protein conformations from M2 calculations. For example, M2 calculations generated 482 distinct conformations of free P1 within 10 RT of the most stable free conformation. Even if free P1 were equally stable in all 482 energy wells with only one bound conformation, the maximum change in conformational entropy would only be reduced by RTln 482 =~3.7 kcal/mol, which is significantly smaller than the -TΔS values in Table 3. We may approximate vibrational entropy through -TΔS vib = -TΔS config + TΔS conf . S3 Table shows that C1 has much smaller vibrational entropy loss than peptides P1-P14. In sum, both conformational entropy and vibrational entropy are attributed to the smaller configuration entropy loss of C1. Interestingly, P7, with a small residue alanine in the P+1 position, has the second smallest entropic penalty in M2 results but not P13, which has two proline residues in the middle of the peptide. P13 managed to partially eliminate the folded conformations because of the geometric constraint proline residues; however, the entropy cost does not decrease substantially due to the big vibrational entropy loss (S7 Fig and S3 Table). Moreover, the restraint by the two prolines resulted in the incorrect orientation of ligand-bound conformations, which significantly weakens the polar attractions (S11 Fig). Inhibitor design: New strategy for promiscuous modular domains?
Two strategies are commonly used in ligand design for enhancing binding affinities: increasing intermolecular attractions and decreasing entropy loss upon binding. For example, new interactions between ligands and receptors, such as adding hydrogen bonds, can be introduced to increase enthalpic attractions [87][88][89][90][91][92][93]. The other way is via reducing the entropy cost by pre-rigidifying the ligand to its bound conformation [94,95]. This pre-organization of the ligand to its bound conformation lessens the decrease in number of rotameric states, and thus affinity is increased primarily because of optimizing the entropic term.
Because the number of potential hydrogen bonds may already be maximized by the presence of the phosphate group, we used the latter strategy to pre-organize a ligand by introducing a benzene ring in the ligand backbone to limit its conformational flexibility. Having a benzene ring in the middle at a certain level prevents the ligand from bending and forming intra-molecular hydrogen bonds like other tetrapeptides do. A new ligand, N1, was synthesized (Fig 2) and its binding to BRCT was tested experimentally. Although the conformations were constrained to some degree to reduce conformational entropy penalty (S7 Fig), the loss from the vibrational part was not reduced enough. The conformational constraints by the benzene ring restricted the ligand rearrangement to optimize the polar and non-polar contacts to the protein, thereby resulting in weak binding (Table 3 and S12 Fig). N1 performed similar to P13, so over-rigidifying a ligand is not advantageous, which suggests the challenge in retaining optimized intermolecular interactions in pre-rigidifying a peptidomimic compound. Previous work in design of potent Cbl(TKB)-binding peptides drew the same conclusion [96]. Therefore, because of conformational flexibility at the binding interface of a modular domain, flexible ligands may be favorable.
Another strategy to lower entropy penalty, although less common, is by introducing a less rigid complex while the molecules bind. Because the strategies to further modify the short peptides to increase their bound conformations may be exhausted, compounds with phosphate groups are a better alternative. On the basis of our calculations and the structure of compound C1, we further modified it to enhance intermolecular attractions by the formation of additional hydrogen bonds between the ligand and BRCT. In the meantime, we kept the template structure intact to maintain its flexibility. We added one hydroxyl group to the para site of the benzene ring of C1, which can form hydrogen bonds with K1690 or N1774 with different bound conformations (S13 Fig). Therefore, designed compound D1 shows improved binding affinity, by 2 kcal/mol, with more negative Δ(U+W) as compared with C1 (Fig 9 and Table 3). As compared with C1, D1 has a stronger Coulombic interaction with BRCT because of the additional hydogen bonds (Table 3). Moreover, because D1 can also adopt mutiple distinct bound conformations, the entropy cost is minimal, as is found in C1. The enthalpy-entropy compensation plot shown in Fig 9 clearly indicates that D1 outperforms other peptides by both increasing intermolecular attraction and reducing entropic penalty.
In summary, designing a pre-rigidified ligand to reduce entropy cost can be tricky considering the potential loss of intermolecular attraction due to lack of proper rearrangement in the bound state. Fortunately, making a ligand more flexible and able to retain its plasticity in the bound conformation provides an effective strategy to reduce entropy cost, while the optimization of interactions between such a flexible ligand and a target protein can further improve binding affinity. Although for designing tight binders such as many drug-protein binding systems, pre-rigidified may still be the best strategy, our study points out a new direction for designing inhibitors targeting promiscuous modular domains and PPIs.

Supporting Information available
Supplemental figures, data, results and detailed experimental method are provided in Supporting Information. All files, including the input and output files for MD simulations, post-analysis and M2 methods, MD trajectories, and results from various energy calculations are freely available upon request (email: chiaenc@ucr.edu).
Supporting Information S1 Table. Sources of initial bound conformations of ligands for MD simulation. Table. Numbers of complex, free ligand and protein conformations from M2 calculation. M2 uses a rigorous conformational search through dihedral distortion for new conformations. Molecular torsional modes are calculated via diagonalization of matrix of energy 2ndderivatives transformed into internal coordinates with all bond and angle rows and columns removed. After a complete distortion along these modes, the whole system is energy minimized via a quasi-Newton geometry optimization to get new conformations. 1RT means numbers of conformations within 1RT above the global energy minimum. (DOCX) S3 Table. Approximated conformational and vibrational entropy (kcal/mol) for P1-P14, C1, N1 and D1. The conformational entropy penalty is approximated through RTln M (M is the number of conformations within 10RT of most stable free ligand conformation from S2 Table). The vibrational entropy penalty was computed by -TΔS vib = -TΔS config + TΔS conf . a For C1 and D1, they have at least three distinct bound conformations (Figs 5 and S13), so the conformational entropy penalty of C1 and D1 is approximated through RTln (M/3).