Probing Difference in Binding Modes of Inhibitors to MDMX by Molecular Dynamics Simulations and Different Free Energy Methods

The p53-MDMX interaction has attracted extensive attention of anti-cancer drug development in recent years. This current work adopted molecular dynamics (MD) simulations and cross-correlation analysis to investigate conformation changes of MDMX caused by inhibitor bindings. The obtained information indicates that the binding cleft of MDMX undergoes a large conformational change and the dynamic behavior of residues obviously change by the presence of different structural inhibitors. Two different methods of binding free energy predictions were employed to carry out a comparable insight into binding mechanisms of four inhibitors PMI, pDI, WK23 and WW8 to MDMX. The data show that the main factor controlling the inhibitor bindings to MDMX arises from van der Waals interactions. The binding free energies were further divided into contribution of each residue and the derived information gives a conclusion that the hydrophobic interactions, such as CH-CH, CH-π and π-π interactions, are responsible for the inhibitor associations with MDMX.


Introduction
Recently, due to key roles in maintaining the integrity of the genome, the tumor suppressor protein p53 has been paid an extensive attention [1]. P53 protein can coordinate the cellular response to DNA damage by inducing cell cycle arrest or apoptosis [2]. Active p53 can efficiently hold back tumor's growth and protect human health [3]. The previous studies suggested that the inactivation of p53 is tightly related with human cancer, which mainly arises from either point mutations in TP53 gene or functional inhibition by negative regulators MDM2 and MDMX [4][5][6][7][8]. The data from clinical treatments proved that over-expressions of two proteins MDM2 and MDMX were found in~50% of all cancer patients around the globe, which significantly influences the wild-type function of p53 [9][10][11].
MDMX is also named as MDM4 and highly homologous to MDM2. It is another significant negative regulator of p53 [4,12]. Three key residues (Phe19´, Trp23´and Leu26´) can form direct interactions with MDM2/MDMX [13][14][15]. Although MDMX has overall similar structure to MDM2, it does not possess the capability of ubiquitin ligase and its expression level is not related with p53 [16][17][18][19]. Several previous experimental studies also demonstrated that inhibition effect of current inhibitors on MDMX is weaker than that on MDM2 [7,[20][21][22][23]. Thus, it is essential to further reveal the interaction mechanism of inhibitors with MDMX for development of potent inhibitors rescuing the wild-type function of p53.
Up to now, many researches have focused on development of inhibitors restoring the p53 function [24][25][26][27][28][29]. Lots of existing inhibitors displayed strong inhibition ability on the p53-MDM2 interaction, but they can not efficiently control the association with MDMX. Detailed clarification of the binding modes of peptide and non-peptide inhibitors to MDMX can contribute valuable information on the structure-affinity relationship of the MDMX binding compound, which is helpful for designs of efficient inhibitors. Joseph et al. applied MD simulations and binding free energy calculations to reveal the cause of the binding difference of p53 and nutlin to MDM2 and MDMX [30]. Li et al. explored the inhibition mechanism of inhibitors on the p53-MDM2/MDMX interaction by performing systematic mutational analysis on p53 and peptide inhibitor PMI [31]. Chen et al. combined MD simulations and computational alanine scanning method to successfully investigate the difference in binding modes of inhibitors to MDM2 and MDMX [32]. Due to high flexibility of hydrophobic cleft of MDMX, further insights into the binding mode and conformational changes of MDMX induced by inhibitor bindings is of importance for the designs of potent inhibitors interrupting the p53-MDMX interaction.
Two peptide inhibitors PMI and pDI, together with two non-peptide inhibitors WW8 and WK23, were picked out to study their binding modes to MDMX. PMI is a peptide inhibitor (TSFAEYWNLLSP) designed by Pazgier et al. This inhibitor not only structurally shares three common residues (Phe3´, Trp7´and Leu10´) with p53, but also has a strong inhibition ability scaled by K d value of 4.15 nM [33]. pDI is another peptide inhibitor (LTFEHYWAQLTS) contributed by Phan et al., which shows a weak inhibition potency on MDMX with IC50 value of 550 uM and structurally also shares three key residues with p53 [34]. WW8 and WK23 are two non-peptide inhibitors studied by Popowicz et al. and their Ki values are 11 and 36 uM [14], respectively. The structures of these inhibitors were displayed in Fig 1. Although pDI and PMI share three same key residues with p53 and the structures of WW8 and WK23 are also highly similar, their binding affinities to MDMX are greatly different. Thus it is of significance to understand the cause leading to these difference of binding abilities for designs of small molecule inhibitors blocking the p53-MDMX interaction.

Starting structures
The crystal structures of three compounds with pDI, PMI and WW8 taken from the protein data bank were used as the starting model of MD simulations and their corresponding PDB entries are 3EQY, 3JZO and 3LBJ, respectively [14,33,34]. Because the crystal structure of the WK23-MDMX compound is not available, it can be got by modifying the WW8-MDMX structure. The Leap program included in Amber was employed to connect all missing hydrogen atoms with the heavy atoms of MDMX [57]. During the initialization of the simulated structures, all water molecules originating from the crystal structures were retained. The force field parameters of proteins and crystal water molecules were gained by using the FF03 force field [58]. The structures of WK23 and WW8 inhibitors were minimized at semi-empirical AM1 level and BCC charges were assigned to them by using the Antechamber program in Amber package. The force field parameters of WK23 and WW8 were generated by using the general Amber force field (GAFF) [59]. To neutralize the charge of each system, an appropriate number of chloride counterions were added to four inhibitor-MDMX compounds. To reflect the solvent environment of proteins, an octahedral periodic box consisting of TIP3P water molecules with a buffer of 12.0 Å was adopted to solvate the initialized system.

MD simulations
To relieve bad conformation between atoms, energy minimization and MD simulations were performed on all initialized systems by using the Amber dynamics program. Firstly, the starting models were minimized in 3000 steps with constraint on the compound using a harmonic restriction of a strength of 100 kcal·mol -1 ·Å -2 . Next, full minimization of another 3000 ps without any restraints were run [32,60]. Then, the temperature of system was slowly heated from 0 to 300 K in 500 ps. To achieve a reliable stability, the system endured an equilibration of another 500 ps at 300K. Finally, 60-ns MD simulations without restriction were conducted on each system at 1 atm and 300 K to obtain a stable trajectory for the following post-processing analysis. The covalent bonds connecting with hydrogen atoms were restrained by using the SHAKE method so that the time step of 2 fs was used during MD simulation [61]. The Particle Mesh Ewald (PME) method was adopted to compute the long-range electrostatic interaction [62,63]. The cutoff distances for the non-bonded interactions were set to 10 Å.

Solvated interaction energy method
The SIE method developed by Naim et al. was usually applied to predict the binding free energies of inhibitors to proteins in many life systems by using the following equation [64]: where E c is the intermolecular Coulomb interactions in the bound state and the parameter D in is the solute interior dielectric constant. The term E vdW represents the van der Waals interaction energies of inhibitors with proteins. These two terms can be computed by using the Amber molecular force field FF03 [58]. The term ΔG R reflects the reaction field energy difference generated by inhibitor associations, and was obtained by solving the Poisson equation using the boundary element method BRI BEM [65,66] and a variable-radius solvent probe [67]. The fourth term γ • ΔMSA(ρ) is related with the change of the molecular surface area induced by inhibitor associations and the parameter γ is the molecular surface area coefficient. The parameters α and C are the global proportionality coefficient related to the loss of conformational entropy upon binding and the constant, respectively. All parameters involved in the current calculations are set to α = 0.1048, D in = 2.25, D in = 2.25, ρ = 1.1, γ = 0.0129 kcal/(molÁÅ2) and C = −2.89 kcalÁmol −1 , respectively [64,68,69]. Although the coefficients for SIE were parameterized using the AMBER ff99SB force field [70], our previous studies on the interaction mechanism of inhibitors with MDM2 based on FF03 force field proved that these coefficients of the SIE method can also produce rational results in binding free energy predictions [71]. The SIE method has obtained great success in insight into the structure-affinity relationship of several inhibitor-protein binding systems [69,72]. Thus, the SIE method was adopted to evaluate the binding abilities of inhibitors to MDMX by using the program Sietraj and 200 snapshots taken from the last 20 ns of MD trajectory at an interval of 100 ps [68].

MM-GBSA calculations
MM-GBSA method [73] can be also used to predict the binding free energies of inhibitors to proteins. The previous studies suggested that this method has obtained success in studying the inhibition mechanisms of the p53-MDM2/MDMX association [32,47,53]. Thus, MM-GBSA method were implemented to compute the binding free energies of the current four inhibitors to MDMX based on the following empirical equation.
in this equation, the first two terms ΔE ele and ΔE vdw represent the electrostatic and van der Waals interactions, respectively and were calculated by using the Amber molecular mechanics force field FF03. The term ΔG pol is the polar contribution to solvation free energy, which was computed by using the modified GB model developed by Onufriev et al. [74]. The fourth term (ΔG nonpol ) in Eq 2 is the non-polar solvation energy and was computed by using the empirical Eq 3 in the current calculation, the parameters γ and β were assigned as 0.0072 kcalÁmol −1 ÁÅ −2 and 0.0 kcalÁmol −1 , respectively [75]. The last term (TΔS) is the entropy contribution and comes from the changes in the translational, rotational and vibrational degrees of freedom induced by inhibitor bindings, which is calculated by using the normal-mode analysis [76]. In this work, 200 snapshots taken from the last 20 ns of MD trajectory at an interval of 100 ps were used for MM-GBSA calculation and 50 snapshots from the last 20 ns of MD simulations at an interval of 400 ps for the normal-mode analysis.

Internal dynamics
The cross-correlation analysis is a significant tools to study the conformation change and internal dynamics of proteins. To probe the changes of internal dynamics of MDMX induced by inhibitor associations, the cross-correlation matrix C ij for C α atoms i and j can be constructed by using the following equation [77].
in which Δr i describes the displacement of the ith atom relative to its average position. The angle bracket indicates an ensemble average over the frames extracted from the equilibrated MD trajectory. The values of C ij fluctuate in the range from -1 to 1. The positive value of C ij indicates a correlated motion between C α atoms, while the negative one implies an anti-correlated motion.

Results and Discussions Equilibrium of dynamics simulations
The root-mean-square deviations (RMSD) of backbone atoms in MDMX relative to the initial minimized structures through the phase of MD simulations were calculated and depicted in Fig 2. These RMSD values can be used to evaluate the reliability of MD simulation equilibrium.   2. indicates that the RMSD values of the compounds with pDI, WK23 and WW8 are higher than that of the PMI-MDMX compound, which implies that the restriction of these three inhibitor bindings on MDMX may be weaker than that of PMI.

Dynamics of MDMX
Root mean square fluctuation (RMSF) is generally applied to evaluate the fluctuation of a specific residues relative to a reference structure. To reveal the flexibility of MDMX, the averaged RMSF values of C α atoms were computed by using the trajectories after the MD equilibrium (Fig 3). To further discover the impact of inhibitor associations on internal dynamics of MDMX, the cross-correlation matrices describing the fluctuations of C α atoms were computed (Fig 4).  Additional anti-correlated motions also appear in the diagonal region consisting of the residues 90-107. Except for the above anti-correlated motions, slight correlated motions can be found in the diagonal region indicated by yellow.
The comparison with the PMI-MDMX compound shows that the associations of pDI, WK23 and WW8 with MDMX lead to obvious changes of internal motions (Fig 4B-4D). The presence of these three inhibitors strengthen not only the correlated motions in the domain (the residues 49-65), but also the correlated motions of the inter-domain between the residues 88-93 and 71-76. For the WW8-MDMX compound, the association of WW8 induces the increase of correlated motions in the diagonal region consisting of the residues 94-107 (red). The bindings of inhibitors WK23 and WW8 to MDMX result in the disappearance of anti-correlated motions of the inter-domain between the residues 65-72 and 51-58, while the pDI binding increases this anti-correlated motions. Finally, the anti-correlated motions between the residues 75-88 and 27-40 disappear due to the restriction of these three inhibitors. The above results prove that different inhibitor bindings induce the changes of internal dynamics in MDMX during MD simulations, which reflects that the flexibility of MDMX is large and inhibitor bindings produce the different dynamic behavior in MDMX. The analysis reported here basically agrees with the previous RMSF analysis.
To evaluate the overall effect of inhibitor associations on the flexibility of MDMX, the eigenvalues were obtained by diagonalizing the covariance matrix constructed by the atomic coordinate [78,79] and the information was displayed in Fig 5. The first several eigenvalues indicate the concerted motions of residues. It is observed that these eigenvalues quickly reduce in amplitude to reach a number of constrained and more localized movement. The first six principal components contribute up to 58.2%, 59.8%, 61.4% and 69.4% in the total motions for the PMI-MDMX, pDI-MDMX, WK23-MDMX and WW8-MDMX compounds, respectively. This result suggests that the bindings of different structural inhibitors produce different effects on the concerted motions of MDMX. The above analyses basically agree with the previous studies [29,30,60].

Binding free energy calculations
To clarify the binding mechanism of inhibitor bindings to MDMX, binding free energies were computed by using MM-GBSA method based on MD trajectories. The computed results were given in Table 1. One can see from Table 1 that the binding free energies of PMI, pDI, WK23 and WW8 to MDMX are -13.3, -10.7, -8.4 and -8.9 kcalÁmol −1 , respectively. It is encouraging that the rank of our predicted binding free energies agrees well with the experimentally determined one. As shown in Table 1, although the predicted binding free energies are about -2.2 kcalÁmol −1 higher than the corresponding experimental values, a good correlation was observed between our predicted values and the experimental ones. Although Cheng et al also computed the binding free energies of PMI and WW8 to MDMX, their predicted results is  Table 1 shows that the van der Waals interaction (ΔE vdw ) and non-polar solvation energy (ΔG nonpol ) favor the association of inhibitors with MDMX. Although the electrostatic interactions (ΔE ele ) provides a favorable force to the inhibitor bindings, it completely screened by the polar solvation energy to lose its advantage. Furthermore, it is noted that van der Waals interaction is much stronger than non-polar solvation energy. Thus, a conclusion was obtained that the van der Waals interactions mainly control the bindings of these inhibitors to MDMX.
To complement our results predicted by MM-GBSA method, the SIE method was also employed to compute the association affinity of these inhibitors to MDMX (Table 2). Table 2 suggests that the binding free energies of PMI, pDI, WK23 and WW8 to MDMX are -10.71, -9.92, -6.89 and -7.09 kcalÁmol −1 , respectively. Moreover, the rank of the predicted binding free energies is also in good accordance with the experimental determined one. The above results from two methods show that the current free energy analysis is reliable.
By comparing the results from two methods, it is seen that the ranks of binding free energies predicted by two methods agree well with the experimental one, but the binding free energies predicted by the SIE method are quantificationally closer to experimental values than those by MM-GBSA method. Based on these two calculations, a common conclusion was derived that the van der Waals interactions are responsible for the bindings of inhibitors to MDMX. By comparison of free energy components, it is found that two peptide inhibitors PMI and pDI can form more interaction contacts than two non-peptide inhibitors. The above analyses basically agree with the previous studies [30,32,60].

The structure-affinity relationship
The inhibitor-residue interactions were computed by using the residue-based energy decomposition method [75] to evaluate the contributions of separate residues to the inhibitor association. The obtained detailed inhibitor-protein interaction spectra (Fig 6), analyses of hydrogen bonding interactions (Table 3) and the positions of the inhibitors relative to the key residues in the MDMX (Fig 7) were combined to probe the structure-affinity relationship of the binding compound. Fig 6A shows that the binding energies of nine residues Thr48, Met53, Ile60, Met61, Tyr66, His72, Val92, Tyr99 and Arg103 from MDMX to PMI are stronger than -1.30 kcalÁmol −1 . The interaction of Thr48 with PMI is about -1.31 kcalÁmol −1 , which may arise from the hydrophobic interaction between the CH groups of Thr48 and the alkyl of Leu10 0 (Fig 7A). Among nine residues, the interaction energy of Met53 with PMI is the strongest (-4.87 kcalÁmol −1 ). Analyses show that this interaction may arise from two contributions: (1) the CH-π interaction between the indole of Trp7 0 and the alkyl of Met53, (2) the hydrogen bonding interaction between the The signs "±" represent standard errors.
doi:10.1371/journal.pone.0141409.t002  (Table 2), with the distance of 2.83 Å and the occupancy of 92.31% between the corresponding oxygen and nitrogen atoms (Table 3). Structurally, the alkyls of Ile60 and Met61 are near the phenyl of Phe3 0 and apt to generate the CH-π interactions, and their interaction intensities with PMI are -1.41 and -1.89 kcalÁmol −1 , respectively. The binding intensity of Tyr66 to PMI is -2.49 kcalÁmol −1 , which structurally agrees with a paralleled π-π interaction between the phenol of Tyr67 and the phenyl of Phe3 0 . The imidazole of His72 is in the vicinity of the phenyl of Phe3 0 to generate a π-π interaction, which provides an energy contribution of -1.65 kcalÁmol −1 . The interaction intensity of Val92 with PMI is -2.81 kcalÁmol −1 , which mainly originates from the CH-π interaction between the alkyls of Val92 and the indole of Trp7 0 (Fig 7A). According to Fig 6A, the residue Tyr99 also forms a very strong  interaction with PMI (-4.47 kcalÁmol −1 ), which is in good agreement with the CH-π interaction of the alkyl in Leu10 0 with the phenol of Tyr99. Additionally, the positively charged residue Arg103 can generate strong polar interactions with the polar residues in the C terminus of PMI, which provides a contribution of -1.93 kcalÁmol −1 to the PMI association. By Comparison with the PMI-MDMX compound, it is observed from Figs 6B and 7B that eight residues in MDMX form strong interactions with pDI. These residues involve Lys50, Met53, Ile60, Met61, Tyr66, His72 Val92 and Tyr99. Except for Lys50, the interactions of the other seven residues with pDI are similar to those with PMI. The interaction energy of Lys50 with pDI is about -4.80 kcalÁmol −1 , which is mainly from the CH-CH interaction of CH group in Lys50 with Leu10 0 . Differently, the interactions of Met53, His72, Tyr99 and Arg103 with pDI decrease obviously, and the interaction of Thr48 with pDI disappears.
Figs 6C and 7C indicate that eight residues in MDMX produce strong interactions with WK23, and these residues are Met53, His54, Leu56, Gly57, Ile61, His72, Val92 and Tyr99. Compared to Fig 6A, the interaction modes of Met53, Ile61, His72, Val92 and Tyr99 with WK23 are similar to those with PMI and only interaction intensity changes. The interaction energy of His54 with WK23 is -1.34 kcalÁmol −1 , which mainly stems from the π-π interaction between the imidazole of His54 and the hydrophobic ring R2 of WK23. The alkyls of Leu56 and the CH group of Gly57 are near the ring R2 of WK23 to contribute the CH-π interaction energies of -1.55 and -2.2 kcalÁmol −1 for the binding of WK23 to MDMX, respectively.
As seen in Figs 6D and 7D, the residues Lys50, Met53, His54 and His95 in MDMX produce obvious interactions with WW8. The interactions of Met53 and Gly57 with WW8 are similar to those with WK23. Structurally, the alkyls of Lys50 is close to the hydrophobic ring R4 of WW8 to generate a CH-π interaction, which contributes an energy of -1.38 kcalÁmol −1 . The energy of WW8 with His95 is -2.31 kcalÁmol −1 , which structurally agrees well with the π-π interaction of the hydrophobic ring in His95 with the rings R1 and R2 of WW8. The above results basically agree the previous studies [14,30,32].
Based on the above analyses, the interactions of separate residues in MDMX with pDI, WK23 and WW8 obviously change relative to those with PMI. The interactions of Tyr67, His72 and Tyr99 in MDMX with PMI are stronger than those with pDI, WK23 and WW8. These changes in the inhibitor-residue interactions correspondingly induce the changes of internal dynamics in MDMX. The previous cross-correlation analysis and RMSF values agree with the results reported here. In summary, the hydrophobic CH-CH, CH-π and π-π interactions drive the bindings of inhibitors to MDMX, thus optimization of these hydrophobic interactions is very helpful for successful designs of potent inhibitors targeting the p53-MDM2 interactions.

Conclusion
60 ns molecular dynamics simulations were performed on four inhibitor-MDMX compounds to probe the effects of inhibitor bindings on the conformational changes of MDMX. The crosscorrelation analysis and RMSF calculation based on MD trajectory show that the flexibility of MDMX binding cleft is large and the dynamic behavior of residues are different due to the presence of different structural inhibitors. The calculations of binding free energies by MM-PBSA and SIE method prove that van der Waals interaction drives the inhibitors binding to MDMX. The inhibitor-residue interactions were computed by using the residue-based energy decomposition method to obtain a detailed insight into the inhibitor-MDMX binding modes. The results indicate that the CH-CH, CH-π and π-π interactions between inhibitors and separate residues are the main forces controlling the inhibitor association to MDMX. This study can provide significant information on the structure-affinity relationship of the binding compound for the development of effective inhibitors targeting the p53-MDMX interaction.