Revealing the binding modes and the unbinding of 14-3-3σ proteins and inhibitors by computational methods

The 14-3-3σ proteins are a family of ubiquitous conserved eukaryotic regulatory molecules involved in the regulation of mitogenic signal transduction, apoptotic cell death, and cell cycle control. A lot of small-molecule inhibitors have been identified for 14-3-3 protein-protein interactions (PPIs). In this work, we carried out molecular dynamics (MD) simulations combined with molecular mechanics generalized Born surface area (MM-GBSA) method to study the binding mechanism between a 14-3-3σ protein and its eight inhibitors. The ranking order of our calculated binding free energies is in agreement with the experimental results. We found that the binding free energies are mainly from interactions between the phosphate group of the inhibitors and the hydrophilic residues. To improve the binding free energy of Rx group, we designed the inhibitor R9 with group R9 = 4-hydroxypheny. However, we also found that the binding free energy of inhibitor R9 is smaller than that of inhibitor R1. By further using the steer molecular dynamics (SMD) simulations, we identified a new hydrogen bond between the inhibitor R8 and residue Arg64 in the pulling paths. The information obtained from this study may be valuable for future rational design of novel inhibitors, and provide better structural understanding of inhibitor binding to 14-3-3σ proteins.

Scientific RepoRts | 5:16481 | DOI: 10.1038/srep16481 Each 14-3-3 proteins consists of characteristic cup-like shape functional dimers with each monomer has nine antiparallel α -helices displaying a so-called amphipathic groove that accommodates the mostly phosphorylated interaction motifs of their partner proteins (see Fig. 1A) 11,12 . Small-molecule regulation on PPIs is one of the most exciting but also difficult fields in drug development and chemical biology 13 .
Previously, several attempts have been made to develop small-molecule inhibitors for the 14-3-3 PPIs. For example, Wu et al. designed and synthesized a peptide-small-molecule hybrid library based on the original optimal 14-3-3 binding peptide and maintained the central phosphoserine residue 14,15 . Corradi et al. employed an in silico structure-based inhibitor design approach to identify the first non-peptidic small molecule compounds with anti-proliferative activity 16 . Zhao et al. identified and experimentally confirmed a pyridoxal-phosphate derivative, which create a covalent linkage of the pyridoxal-phosphate moiety to the residue Lys120 in the binding groove of the 14-3-3 protein 17,18 . Bier et al. reported a molecular tweezers which bind to a 14-3-3 adapter protein and modulate its interaction with the partner proteins 19 . Thiel et al. identified noncovalent and non-peptideic small-molecule inhibitors for extracellular 14-3-3 PPIs by virtual screening 20 .
In the work by Thiel et al., the crystallographic structures of the 14-3-3σ protein and inhibitors complexes were solved. Such high-quality structural data can be exploited to design the PPI inhibitors in silico 20 . This is very important for the understanding the protein-inhibitor interactions at the atomic level of this class of compounds, which may lead to the development of 14-3-3σ inhibitors with better potency. It is well-known that molecular dynamics (MD) simulations can enhance our understanding of binding mechanisms for protein-inhibitor complexes, such as the 14-3-3σ protein and its inhibitors complexes, by providing quantitative binding affinities [21][22][23][24][25][26][27] . Several computational methods with various levels of computational expense and accuracy can be used to estimate the inhibitor binding affinities and selectivities. These methods include the thermodynamic integration (TI), the free energy perturbation (FEP) method 28,29 , and molecular mechanics generalized Born surface area (MM-GBSA) method 30,31 . Among them, MM-GBSA method is a versatile tool for calculating the binding free energy of a given protein-inhibitor complex. In this method, the gas-phase energy, calculated using conventional molecular mechanics force fields such as AMBER 32 , is combined with a continuum model of solvation that includes a surface area based nonpolar contribution 33 and a polar solvation free energy calculated with the generalized Born (GB) approximate model of electrostatics 34 . Noted that MM-GBSA method utilizes a fully pairwise potential to decompose the total binding free energy into atomic/group contributions in a structurally nonperturbing formalism 30 .
Steered molecular dynamics (SMD) simulation takes inspiration from single-molecule pulling experiments 35 , and dissociates a complex structure by a pulling force 36,37 . The non-equilibrium dynamics of the system under a pulling force can map out the free-energy landscape in terms of the potential of mean force (PMF) 38 with high precision and efficiency [39][40][41][42] . The free-energy difference between the bound states and the dissociate states can be extracted by measuring the work along the transition paths. Thus, SMD simulations have become widely used in studying biochemical processes including the unfoulding/ foulding mechanism of proteins 43 , transportation of ions and organic molecules across membrane channels 39,[44][45][46][47] , and the mechanisms of protein-inhibitor binding 36,40,48 .
In this paper, we combined the MD simulation with MM-GBSA method to calculate the binding free energies between the 14-3-3σ and its eight inhibitors (Fig. 1B). Our calculated binding free energies are in agreement with the experimental results. The interaction between the phosphate group of inhibitors and the hydrophilic residues are the main contribution for the binding free energies in all compounds (14-3-3σ proteins and inhibitors). To explore the unbinding mechanism for 14-3-3σ and its inhibitors, SMD simulations combined with Brownian-dynamics fluctuation-dissipation theorem (BD-FDT) were used to calculate the interaction energy landscape of 14-3-3σ with inhibitors R1 and R8. Base on the binding model of 14-3-3σ and its inhibitors, a new inhibitor R9, which can form a new hydrogen bond between group R9 = 4-hydroxypheny and residue Glu57, was designed and evaluated in this work.

Results
The protonation of the phosphate group. The crystallographic complex of the phosphate peptide and the 14-3-3σ protein (PDB ID: 1YWT) 49,50 revealed that the phosphate group of the binding peptide forms several hydrogen bonds with 14-3-3σ protein. The structure-based net charges at neutral pH for the 14-3-3σ protein were calculated by using the Adaptive Poisson-Boltzmann Solver (APBS) and PDB2PQ program 51 and visualized resulting electrostatic potentials in VMD software 52 (Fig. 2A). It is clear that the groove in 14-3-3σ protein is hydrophilic 53 . The hydrophilic pocket of the phosphate group is formed by several hydrophilic residues (Arg60, Arg133, Tyr134 and so on). In our previous work, the phosphate group in phosphoserine residue was in unprotonated state 50 . So we set the phosphate group of Stability of the compounds. MD simulations for eight compounds were performed for the time duration of 20 ns. The root mean square deviations (RMSDs) from the crystallographic structure, which can effectively assess the dynamics stability of compounds, were analyzed by using Ptraj 54 module of AmberTools software for apo-14-3-3σ , as well as for compounds R1 and R8 (see Fig. 3). The average RMSDs of binding pocket in the last 5 ns MD simulations for apo-14-3-3σ (1.62 ± 0.16 Å) is larger than that for compound R1(1.17 ± 0.14 Å), as well as that for compound R8 (1.15 ± 0.13 Å). This indicates that the binding pocket of the 14-3-3σ protein is more stable with inhibitor than that without inhibitor. It is noted that the RMSDs for the inhibitors show large fluctuation (Fig. 3), indicating some groups of the inhibitor would not bound tightly to the proteins. To evaluate which part of the inhibitor fluctuate largely, we extracted two groups (group one: 2-hydroxyphenylphosphonic acid; and group Rxs: which names are shown in Fig. 1B) of inhibitors to calculate their RMSDs. The standard deviations of the RMSD for the inhibitors (0.33 Å and 0.55 Å) are larger than those for the group one (0.18 Å and 0.27 Å) and smaller than those for the group Rxs (0.52 Å and 0.74 Å) for compounds R1 and R8, respectively, as well as for compounds R2-R7.
Analysis of binding free energy. We noted that 14-3-3σ would undergo conformational change caused by the binding to inhibitor. However, since we are more concerned with the ranking of the calculated binding free energies for all inhibitors with the same chemical scaffold (Fig. 1B), all the snapshots used in the MM-GBSA were extracted from the trajectories of the compounds. The binding free energies for all eight systems were calculated by using mm_pbsa program in AMBER 12 and summarized in Table 1. Though the predicted absolute free energies were larger than those of the experimental results, the ranking orders of them were in good agreement. Figure 4 shows how well the predicted free energies reproduce the experimental data. The correlation coefficient r 2 is 0.93. Besides ranking order of the binding free energies correctly, MM-GBSA method can decompose the total binding free energy into individual components, thereby enabling us to understand the complex binding process in detail 31 . For the eight compounds, the van der Waals interactions and the nonpolar solvation energies, which are responsible for the burial of inhibitor's hydrophobic groups upon binding, are favorable for binding free Identification of the key residues responsible for the binding of inhibitor. In order to find which residues make significant intermolecular interaction contributions to the binding with the inhibitors, the decomposition of the electrostatic interaction energy, van der Waals energy, and solvation free energy for all compounds were analyzed and the results are depicted in Fig. 5 for compounds R1 and R8 and in Fig. S1 for compounds R2-R7, respectively. The decomposition method with MM-GBSA can naturally be used for the energy decomposition at the atomic level for the per-atom contributions summed over all atoms of each residue to obtain the contribution of each residue. This has been successfully applied to a lot of protein-inhibitor binding systems. The major favorable energy contributions originate predominantly from seven residues (Lys53, Arg60, Lys126, Arg133, Tyr134, Leu178, and Val182) with averaged energy contribution larger than − 0.5 kcal/mol in all compounds. Special attention had been paid to three residues (Arg60, Arg133 and Tyr134) with large electrostatic contribution. For example,  the electrostatic contributions of residues Arg60, Arg133 and Tyr134 are − 17.37, − 19.04, and − 5.0 kcal/ mol for compound R1, respectively. The phosphate group has negative charge and residue arginine has positive charge, resulting in strong electrostatic attraction between them. The hydrogen bonds between the phosphate group and the 14-3-3σ protein were listed in Table 2, showing the occupancies and distances of hydrogen bonds in all compounds. The phosphate group forms three hydrogen bonds with both residues Arg60 and Arg133, as well as one hydrogen bond with residue Tyr134. Most of the hydrogen bonds are stable with high occupancy and similar distance in all compounds ( Table 2), implying that the phosphate groups were tightly bonded in the binding pocked formed by three hydrophilic residues (Arg60, Arg133 and Tyr134). This result is in accordance with the analysis of RMSDs. The side chain of residue Lys53 is in the binding pocket of the phosphate group and may contribute large electrostatic interaction energy. However, the total contributions for Lys53 in compounds R5, R6,   R7, and R8 are less than 1.0 kcal/mol, which are less than those in compounds R1, R2, R3, and R4. Although the gas-phase electrostatic interaction of Lys53 is stronger in compounds R5, R6, R7, and R8, it is compensated by the polar solvation energy. We calculated the averaged distances between the nitrogen atom of side chain of Lys53 and the phosphorus atom of the phosphate group over the last 5 ns MD trajectories. By contrast, the averaged distances are smaller in compounds R1 (3.60 Å), R2 (3.57 Å), R3 (3.85 Å), and R4 (5.21 Å) than those in compounds R5 (7.66 Å), R6 (6.30 Å), R7 (6.13 Å) and R8 (8.11 Å). In order to understand the local structural features between the residue Lys53 and inhibitors in compounds R1 and R8, their relative position are shown in Fig. 6A,B, respectively. It is clearly seen from Fig. 6B that there are a few water molecules between the phosphate group of R8 and the side chain of Lys53. As shown in Fig. 6C, there are strong interactions between three residues (Lys126, Leu178, and Val182) and group one of R1, as the van der Waals energies are favorable the binding for residues Leu178 and Val182 to group one, while the electrostatic energies for Leu126. There is a π -alkyl interaction (− 0.69 kcal/mol) between the side chain of Val182 and the group one of inhibitor R1. There are three unfavorable residues (Asp130, Glu137, and Glu186) for inhibitor binding to protein. The averaged free energies for these three residues in eight compounds are 0.93, 1.03, and 0.97 kcal/mol, respectively. These free energies also attributed to the electrostatic interaction. Since the residues aspartic acid and glutamic acid have negative charges, they repel the phosphate group and attract the residues with positive charge in the binding pocket. It is clear that the key residues mainly interact with the group one of the inhibitors, resulting in the formation of a pocket surrounding the group one ( Fig. 2A). This is in agreement with the state that the phosphate has the strongest pharmacophoric properties 20 . By contrast, the Rxs are surrounded by several residues, while there is no stronger interaction between Rxs and residues ( Fig. 5 and Fig. S1).
SMD simulation combined with BD-FDT. SMD simulations were performed to investigate the dynamic processes of two inhibitors (say R1 and R8) unbinding from the 14-3-3σ protein. The starting structures of compounds R1 and R8 for SMD simulations were extracted from the last structure of the afore-presented MD simulations. Then the starting structures were rotated for the orifice of the inhibitor binding pocket toward the + z direction, put them in a box of water, and neutralized the systems. Then 10 ns equilibrated MD simulation was carried out for each system. In our SMD simulations, each inhibitor is represented by two centers. Both centers were steered at the same time along z direction. The pulling speed was set at 0.01 Å/ps in z direction. In order to reduce the impact of pulling on the 14-3-3σ protein, the inhibitor can move freely in x and y directions, and the whole pulling path was divided into 16 segments along the z-direction with 1 Å for each segment. One pulling path way of compound R1 was show in Fig. 7A, the displacement is 16 Å from the bound state to the dissociated state, as well as 25 Å in the xy plane. For each segment, two types of SMD simulations were performed: one for pulling back to (denoted as reverse) the binding site and one for pulling away (denoted as forward) from the binding site. We sampled four forward and reverse pulling paths during which the work done to the system was recorded for each segment. The curves of works done to the systems along the pulling paths are shown in Fig. S2. From these works, we calculated the PMFs as a function of the displacement of inhibitors along z-axis by using the BD-FDT and the results are shown in Fig. 7B. We can see that the PMF difference between the bound state to the dissociated state are − 13.88 and − 9.24 kcal/mol for compounds R1 and R8, respectively. For compound R1, the PMF rises all the way until the displacement reaches to 7 Å where the inhibitor is steered out of the binding pocket. After that, the PMF reaches a plateau. For compound R8, the PMF rises with the displacement < 3 Å and reaches an interesting intermediate state around the displacement within 3.5 Å to 4.5 Å. After that, the PMF rises again and then levels off after 6 Å, indicating the inhibitor is in the dissociate state.

Discussion
Our analysis based on RMSDs shows that the group one is more stable than the group Rxs in each compound. This result indicated that the group one bound tightly to the 14-3-3σ protein, and the group Rxs may not. Our analysis on the binding free energies shows that the electrostatic contribution plays important role in the binding of the inhibitor and the 14-3-3σ protein because the group one of inhibitors has negative charge and the bottom of the binding pocket is hydrophilic. The hydrophilic residues (Arg60, Arg133 and Tyr134) at the bottom of the binding pocket formed seven stable hydrogen bonds with the phosphate group and contributed large electrostatic energy. Additional, two residues (Leu178 and Val182) at the binding pocket of group one contribute large van der Waals energies and the residue Leu126 with large electrostatic energies. So the group one is stable in these compounds. The entropic contribution also plays important role, this is in accordance with the large RMSDs fluctuation of the group Rxs.
As the afore-presented results, the hydrophilic interactions play important role in the inhibitor binding. So we calculated the number of hydrogen bonds between the inhibitor and the 14-3-3σ protein from our SMD simulations (shown in Fig. 7C). The number of hydrogen bonds decreases gradually along all the way until 6 Å for compound R8, indicating that there is no new hydrogen bond formed during the SMD process. For compound R1, the number of hydrogen bonds levels off during the first 2 Å of the displacement, showing stability of the hydrogen bonds formed in the starting structure of SMD simulation (Fig. 8A). From 2 Å to 4 Å, an obvious decrease indicates that some hydrogen bonds were broken, which have shown three styles for four forward and reverse SMD simulations. The broken hydrogen bonds were found for the phosphate group with either the residue Arg133 (Fig. 8B) or the residue Tyr134 (Fig. 8C), as well as both (Fig. 8D). After that, the number of hydrogen bonds increases obviously, indicating that there is a new hydrogen bond formed, then gradually decreases to zero at 9 Å. As shown in Fig. 8E, the new hydrogen bond was formed between the phosphate group and the residue Arg64. As the inhibitor was pulled away from the binding pocket, the last hydrogen bond formed in the starting structure was also broken (Fig. 8F).
Hydrogen bonds is important contributor to the specificity of receptor inhibitor interactions. As the bottom of the binding pocket around the group Rx are hydrophilic, it would be possible to form hydrogen bond between the group Rx of the inhibitor and the protein (Fig. 2A). The distance between hydrogen atom of the group Rx and the oxygen atom of backbone of residue Gly57 is 2.49 Å (Fig. 6C). So we can design an inhibitor termed as R9 with group R9 = 4-hydroxyphenyl. The same study for compound R9 were done as for compounds R1-R8. A hydrogen bond, which is our objective to form in compound R9, was observed between resiude Gly57 and the hydroxy of 4-hydroxyphenyl with occupancy of 68.4% and averaged distance of 2.92 Å. The free energy decomposition indicated that the contribution of residue Gly57 is − 0.68 kcal/mol, which is the largest value among all compounds. However, the total binding free energy for compound R9 is − 26.01 kcal/mol, which is smaller than that for compound R1 and larger than those for compounds R2-R8. The hydrogen bond analysis shows that most of the hydrogen bonds of the phosphate group in compound R9 are weaker than those in compound R1. The free energy decomposition also demonstrated that the binding free energies of key residues especially for Arg60 and Arg133 are smaller in compound R9 than those in compound R1. The added hydrogen bond in compound R9 increases the interaction energy of group R9 and decreases the interaction energy of group one.
In conclusion, the phosphate group of the inhibitor is in unprotonated state. The ranking order of the calculated free energy by MM-GBSA method is in agreement with the experimental data. We found that the phosphate group of the inhibitor forms strong hydrogen bonds with three residues (Arg60, Arg133 and Tyr134). The group one is more stable than group Rx in each compound. The electrostatic and entropic contributions are important for the 14-3-3σ protein binding to the inhibitors. We designed a new inhibitor R9. Although an additional hydrogen bond (group R9 and Gly57) was found in compound R9, its binding free energy is only ranked the second among all the inhibitors. We performed two pulling experiments, and found that the residue Arg64 is important in the unbinding paths of compound R8, because this residue forms a new hydrogen bond with the inhibitor R8 in the pulling paths. As a result, our study can theoretically provide dynamics information and guidance for the design of new potent inhibitors targeting the 14-3-3σ protein.

Materials and Methods
System setups. The crystal structures of all compounds determined by Christian Ottmann et al. were used as the starting structures in our MD simulations 20 . Missing loops were obtained from the crystal structure of 14-3-3σ (PDB ID: 3MHR) 55 . All crystallographic water molecules were retained in the starting model. The standard AMBER force field (FF03) 56 was used to describe the protein parameters and water molecules. Single-point calculations with Gaussion 03 were performed to obtain the electrostatic potential around each compound by using Hartree-Fock/6-31 G* basis set 57 . Atomic partial charges of the inhibitors were fitted to the electrostatic potential (ESP) in this study by using the RESP method with the Antechamber module of AMBER12 package 58 . All compounds were solvated in a rectangular periodic box of TIP3P 59 water molecules with a margin distance of 12 Å, and the systems were neutralized by adding an appropriate number of sodions.

Molecular dynamics simulation.
For all systems, energy minimizations and MD simulations were carried out using the AMBER12 package 58 . The solvated models were first minimized in constant volume by 1000 cycles of steepest descent minimization followed by 1000 cycles of conjugated gradient minimization. After energy minimization, the harmonic restraints with force constants of 2 kcal/(mol·Å 2 ) were applied to all atoms of compound and constant volume was carried out for 70 ps, during which the systems were heated from 0 K to 300 K. Subsequent a constant-pressure MD was used for 90 ps to adjust the solvent density. Finally, a 20 ns production run of constant-pressure MD simulation was carried out at 300 K without any restraints. The MD simulations were performed with periodic boundary conditions at 300 K with the Langevin thermostat and long-range electrostatic interactions with particle mesh Ewald method 60 . All the covalent bonds involving hydrogen atoms were constrained by applying the SHAKE algorithm 61 . The time step of all MD simulations was set to be 2 fs with a cutoff of 12 Å. The pressure was kept at 1.0 atm using isotropic positional scaling. The intermediate structures were saved at every 1 ps for analysis.
Free energy calculations. The stable MD trajectory obtained for each compound was used to estimate the binding free energies (∆G bind ) by using MM-GBSA method, which has been implemented in AMBER12 program. A total number of 500 snapshots were chosen evenly from the last 5 ns MD trajectories with an interval of 10 ps. All counterions were stripped, as well as all water molecules. Our MM-GBSA calculations for each snapshot were carried out in the same way as in other protein-inhibitor systems. Briefly, the MM-GBSA method can be conceptually summarized by the following equations: where G complex , G protein and G ligand are the free energies of the complex, receptor and inhibitor, respectively. The binding free energy (∆G bind ) is evaluated by a sum of the changes in the molecular mechanical (MM) gas-phase binding energy (∆E MM ), the solvation free energy (∆G solv ) and entropic (− ∆ T S) contribution. ∆E MM is further divided into a van der Waals (∆E vdW ) and a gas-phase electrostatic energies (∆E ele ). These energies were computed using the same parameter set as that used in the MD simulations. And the solvation free energy (∆G solv ) is further divided into a polar (∆G pol ) and a nonpolar (∆G nonpol ) component. The polar component (G pol ) was calculated with a GBSA module of the AMBER12 suite. The nonpolar component (∆G nonpol ) was determined using ∆ = γ + β G S ASA nonpol , where SASA is the solvent-accessible surface area that was determined with MSMS program 62 with a probe radius of 1.4 Å, γ and β were set to be 0.005 kcl·mol −1 ·Å −1 and 0.0 kcal/mol, respectively 63 . The conformational entropy contributions to the binding free energies were estimated for 125 snapshots using the normal-mode analysis with AMBER NMODE module 64 . Steered molecular dynamics (SMD) simulations. The MD simulations and SMD simulations were performed by using NAMD 2.8 package 65 . We employed periodic boundary conditions in all three dimensions and utilized particle mesh Ewald (PME) method for treating long-range electrostatic interactions. The cut-off distance for long-range interactions was set to 12 Å with as witching distance of 10 Å. The time step was 2 fs for short-range and 4 fs for long-range interactions. The covalent bonds of all hydrogen atoms were fixed to their equilibrium length by SHAKE algorithm 61 . The temperature was maintained at 300 K using Langevin thermost at with the damping constant chosen as 5/ps and pressure at 1 bar. In our SMD, a inhibitor was pulled from a bound state to a dissociated state. The center of mass of the inhibitor was initially set at the origin of coordinates. To maintain the location of the 14-3-3σ protein during the SMD simulations, we restrained the z coordinates of several backbone carbon atoms of the 14-3-3σ protein which are far away from the binding pocket. Two compounds were equilibrated for 150000 steps for each segment without pulling forces and for 100000 steps with pulling forces.
Potential of mean force (PMF) from steered molecular dynamics (SMD) simulation. For each segment, the two end states are denoted as A and B, respectively. The Gibbs free energy difference between state A and an intermediate state r (namely, the PMF) was given by the Brownian-dynamics fluctuation-dissipation theorem (BD-FDT) 38 where K B and T are the Boltzman constant and the absolute temperature, respectively. The brackets with subscript F and R represent the statistical average over the forward and reverse paths, respectively. → W r A is the work done to the system along a forward path when the inhibitor is steered from A to r.
B r is the work done to the system along a reverse path when the inhibitor is steered from r to A.